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We study the dynamics of thermalization and the approach to equilibrium in the classical <fi 4 
theory in 1 + 1 spacetime dimensions. At thermal equilibrium we exploit the equivalence between 
the classical canonical averages and transfer matrix quantum traces of the anharmonic oscillator to 
obtain exact results for the temperature dependence of several observables, which provide a set of 
criteria for thermalization. In this context, comparing to the exact results we find that the Hartree 
approximation is remarkably accurate in equilibrium. The non-equilibrium dynamics is studied by 
numerically solving the equations of motion in light -cone coordinates for a broad range of initial 
conditions and energy densities. The long time evolution is described by several distinct stages, all 
characterized by a cascade of energy towards the ultraviolet. After an initial transient stage, the 
spatio-temporal gradient terms become larger than the nonlinear term, and there emerges a stage 
of universal cascade. This cascade starts at a time scale to independent of the initial conditions 
(except for very low energy density). During this stage the power spectra feature universal scaling 
behavior and the front of the cascade k(t) moves to the ultraviolet as a power law k(t) ~ t a with 
a < 0.25 an exponent weakly dependent on the energy density alone. The wake behind the cascade 
is described as a state of Local Thermodynamic Equilibrium (LTE) with all correlations being 
determined by the equilibrium functional form with an effective time dependent temperature T c ff (t), 
which slowly decreases with time as ~ t~ a . Two well separated time scales emerge: while T e g(t) 
varies slowly, the wavectors in the wake with k < k(t) attain LTE on much shorter time scales. This 
universal scaling stage ends when the front of the ultraviolet cascade reaches the cutoff at a time 
scale ti ~ a~ a . Virialization starts to set much earlier than LTE. We find that strict thermalization 
is achieved only for an infinite time scale. 
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I. INTRODUCTION 

The dynamics of thermalization and relaxation in a field theory is a topic of much current interest both in early 
cosmology as well as in ultrarelativistic heavy ion collisions, where the current experimental program at RHIC and 
the forthcoming LHC will search for a new state of matter, the quark gluon plasma. 

Such current interest on the non-equilibrium dynamics of relaxation and thermalization within the settings of 
cosmology and heavy ion collisions in fact rekindles the issue of thermalization in non- linear field theories. Pioneering 
work in this area was initiated by Fermi, Pasta and UlamQ] (FPU) who were the first to address the question of 
ergodicity focusing on the non-linear dynamics in discrete chains of anharmonic oscillators using one of the first 
computers. Since then this problem has been studied within a variety of models^, 13 with the goal of answering 
fundamental questions on ergodicity, equipartition and in general the approach to equilibrium in non-linear theories 
with a large but finite number of degrees of freedom. 

Within both contexts a detailed understanding of the dynamics and the time scales for thermalization are extremely 
important. In cosmology the inflationary paradigm currently being tested by precise observations of the cosmic 
microwave background assumes that after the inflationary stage a period of particle production and relaxation leads 
to a state of local thermal equilibrium thus merging inflation with the standard hot big bang cosmology^Q. 

In ultrarelativistic heavy ion collisions, the current theoretical understanding suggests that after the nucleus-nucleus 
collision, partons (mainly gluons) are liberated and parton-parton scattering leads to a state of thermal equilibrium |(|. 
It has been recently argued0 that there are several different stages that ultimately lead to thermalization after a 
nucleus-nucleus collision. During most of the non-equilibrium evolution the dynamics is determined by the classical 
Yang-Mills equations because the gluon distribution function is non-perturbatively large ~ l/a s up to a saturation 
scale|6j, |7J. This observation motivates a numerical program to study the early stages of ultrarelativistic heavy ion 
collisions in terms of the non-equilibrium dynamics of classical Yang-Mills fields 

Inflationary scenarios lead to particle production either via parametric amplification of fluctuations (in the case of an 
oscillating inflaton) of the inflaton or by spinodal instabilities during phase transitions|9L llCt ll l| . In both cases the non- 
equilibrium dynamics is non-perturbative and results in a large population of soft quanta whose dynamics is nearly 
classical. The non-equilibrium dynamics of particle production and eventual thermalization are non-perturbative 
and the resulting fluctuations contribute to the evolution of the scale factor, namely the backreaction from the 
fluctuations becomes important in the evolution of the cosmological space-time 10] . Both parametric amplification 
or spinodal decomposition lead to non-perturbative particle production of a band of wavevectors, typically for soft 
momenta^ ^J. This non-perturbatively large population allows a classical treatment of the non-equilibrium 
evolution. 

More recently and again motivated by cosmology and heavy ion collisions a fairly intense effort has been devoted 
to studying thermalization and in general relaxation in classical as well as quantum field theories ^2~!I3- Several 
approximate methods have been implemented, from homogeneous and inhomogeneous mean field (Hartree) ll2l to a 
variety of schemes that include quantum correlations after suitable truncations of the hierarchy of eauations|l3l IbH 
IT5I [To . Il7j| and comparisons between them. A numerical study of the non-equilibrium dynamics in the classical 4> 4 
field theory in 1 + 1 dimensions and a comparison of the classical evolution to that obtained from several approximate 
schemes has been reported in refs. [l3HT6l Uj\ . 
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In this article we study the non-equilibrium dynamics leading to thermalization in the classical <p theory in 1 + 1 
dimensions. As mentioned above the initial stages of non-equilibrium dynamics either in cosmology or in ultrarela- 
tivistic heavy ion collisions is mainly classical. Classical field theory must be understood with an ultraviolet cutoff, or 
equivalently an underlying lattice spacing, to avoid the Rayleigh-Jeans catastrophe. Both in cosmology as well as in 
heavy ion collisions there are natural ultraviolet cutoffs: in cosmology it is the Planck scale (although the distribution 
of particles produced during inflation is typically concentrated at much smaller comoving scales in wavenumbers 10J), 
in ultrarelativistic heavy ion collision it is the saturation scale @. 

Our study of the non-equilibrium dynamics in 1 + 1 dimensional classical field theory is inspired by the description 
of the early stages of the non-equilibrium dynamics in ultrarelativistic heavy ion collisions. The study of the non- 
equilibrium dynamics of a classical field theory is interesting and important on its own right and the 1 + 1 dimensional 
case is the simplest scenario. 

The advantage of studying the classical field theory is that the equation of motion can be solved exactly. 

In fact this is the main reason why it has been studied previously as a testing ground for different approximation 
schemes. 

The goals of this article: 

Our motivation, focus and goals in this article are very different from those of previous studies of the classical </> 4 
theory in 1 + 1 dimensions reported in refs. 0,0, 0,0,0]. While refs.0, El 0, 0, 0] studied whether ergod icitv 
and thermalization are achieved and tested a variety of approximate schemes useful in the quantum case fil^ Tl6l Il7| . 
we seek to provide a more detailed understanding of the main dynamical mechanisms that lead to thermalization 
and to study the approach to thermalization by several different observables. We address several questions on the 
equilibrium and non-equilibrium aspects: i) what are the criteria for thermalization in an interacting theory?, ii) what 
is the mechanism that leads to thermalization?, iii) what is the dynamics for different observables and how they reach 
thermal equilibrium?. 

We focus our study on these issues within the context of classical field theory, which is an interesting and timely 
problem all by itself. It is also likely to describe the initial stages of evolution strongly out of equilibrium at least for 
the relevant cases of very large occupation numbers. 

Main results: The main results of this work are the following: 

• We obtain a series of exact results in thermal equilibrium by exploiting the equivalence between the classical 
equilibrium canonical averages and the quantum traces of the transfer matrix of the anharmonic oscillatorj22| . A 
body of results on the spectrum and matrix elements of the quantum-mechanical anharmonic oscillator |24i l25j 
as well as new calculations allows us to extract exact results for high, intermediate and low temperatures for 
the observables. We compare these exact results to the Hartree approximation for high and low temperatures 
and we find a remarkable agreement. These exact results provide a set of stringent criteria for thermalization. 

• We implement a light cone method to study the classical evolution which is very accurate and stable, maintains 
the underlying Lorentz symmetries and conserves energy exactly (that is to machine accuracy on a computer). 

• For a wide range of initial conditions and energy densities we find several distinct stages in the dynamics: the 
first stage corresponds to initial transients with large fluctuations during which energy is transferred from low 
momentum to higher momentum modes. This first stage ends at a time scale £j ~ 200 in dimcnsionlcss units. 
[That is, choosing both the mass and the coupling equal one, as one can always do in classical field theory]. This 
value for tj turns out to be independent of the lattice spacing and the initial conditions provided the energy 
density E/L is not too small: E/L > 10. After a crossover between the and the </>' 2 terms, a second stage 
emerges during which energy transfer from low to high momentum modes becomes very effective and results in 
an ultraviolet cascade, namely the power spectrum of the field and its canonical momentum acquire support 
for ever larger values of the wavevectors. This ultraviolet cascade leads to a very efficient transfer of energy 
between the interaction term <^> 4 , the (ft 2 and the 7r 2 terms which all diminish and the spatial gradient term (f>' 2 
which grows. After a time scale to ~ 50000 a third stage emerges during which the ultraviolet cascade features 
a universal behavior independent of the lattice spacing and the initial conditions. During this third stage the 
transfer of energy to higher momentum modes, namely the cascade, develops a front k(t) which moves towards 
the ultraviolet cutoff and leaves in its wake a state of local thermodynamic equilibrium characterized by 
an effective temperature T c g(t). That is, two well separated time scales emerge : a short time scale leading 
to a state of local thermodynamic equilibrium with an effective temperature that decreases and approaches the 
equilibrium value on a long time scale. While T e s(t) varies slowly with time the fc-modes below the front of the 
cascade k(t) adjust to local thermodynamic equilibrium with T c s(t) on a much shorter time scale. 
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The power spectrum of the canonical momentum features a universal scaling behavior. This universal scaling 
stage ends when the front of the cascade reaches the ultraviolet cutoff at a time scale t\ ~ (L/E) a -1 /" with L 
the size of the system, E the conserved energy and a the lattice spacing, and with a <~ 0.21 — 0.25 a cutoff- 
independent exponent weakly dependent on E/L. During this stage with t\ » t 3> to the front of the 
cascade evolves as k(t) ~ (E/L) 1 / 4 t a and all observables attain their thermal equilibrium values corresponding 
to the effective temperature T e g(t). At the same time T e g(t) decreases as T e g(t) ~ (tt/2) (E/L) 3 / 4 t~ a . For 
t > t\ thermalization continues further, however the front of the cascade is near the cutoff and universal scaling 
behavior does no longer hold. Exact thermalization is only achieved at infinite time when the power spectrum is 
flat and sharply falls off at the cutoff scale and T e fj(oo) = T coincides with the canonical equilibrium temperature. 

• In agreement with most previous studies we find that thermalization does occur but strictly speaking only at 
infinite time. Moreover, an effective local thermal equilibrium description in terms of a (slowly) time 
dependent temperature T e g(t) is available after the second stage at which the interaction term becomes small as 
compared to the spatio-temporal gradient terms. It should be stressed that during the stage of universal cascade 
we obtain a cutoff- independent expression for T e g(t) and therefore for all observables in effective local thermal 
equilibrium. This is therefore the quasi-equilibrium state of the continuum field theory for times beyond to. 

In section II we summarize the classical (f> 4 theory in thermal equilibrium, obtain exact results for low, intermediate 
and high temperature and compare to the Hartree approximation. In this section we provide a set of stringent criteria 
for thermalization and based on the equilibrium results we discuss the main features expected from the dynamics. 
In section III we introduce the light cone approach to study the dynamical evolution in the cutoff theory. In this 
section we discuss a wide range of initial conditions and energy densities along with the averaging procedure and the 
observables to be studied. In section IV we study in detail the power spectra and the energy cascade. Here we discuss 
the stage during which the power spectrum of the canonical momentum features universal scaling behavior, which we 
refer to as the universal cascade. We establish the regime of validity of this universal cascade and show that the wake 
behind the front of the cascade describes a state of local thermodynamic equilibrium at a temperature T e g(t) which 
slowly decreases towards the equilibrium value. 

Our conclusions as well as several comments and discussions on the continuum limit, the quantum theory and the 
kinetic description are presented in section V. 

II. THERMAL EQUILIBRIUM 

We present here the equilibrium statistical mechanics of the classical <j) 4 theory in 1 + 1 dimensions. 

A. Generalities 

The Lagrangian density of the <p 4 field theory reads 

C m ,\ (v) = ^d^tpd^tp -^-tp 2 - ^tp 4 
and leads to the classical equation of motion 

d^p + m 2 p + A ip 3 = . (2.1) 

The theory is defined in a finite volume with periodic boundary conditions. In 1 + 1 space-time dimensions the 
coupling A has dimensions of (mass) 2 . 

In the classical theory one can always rescale variables absorbing the mass m and the coupling A in the field and 
coordinates in order to define the following dimensionless quantities 

771 

t = mx° , x = mx x , p(x° , x 1 ) = —= <fr(x, t) (2-2) 

v A 

so that the Lagrangian takes the form 

An,A W) = ^ A,l i4>) = ^ (i? - - ^ \ ^ (2.3) 
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where <j) = d<f)/dt and <j>' = d<fi/dx. Using the scaled Lagrangian C\.\ (</>) and the canonical momenta II 
tt — d>, we find for the Hamiltonian 



m 



H[n,<f>] = 



dx 



I 



ip and 



(2.4) 



We will consider periodic boundary conditions (PBC) on a ring of (dimensionless) length L, namely 4>(x + L,t) 
<f>(x,t). 

The equations of motion in terms of the dimensionless field is given by 



= 



In terms of Fourier mode amplitudes, 



dx 



-ikx 



/o Vl 

the Hamiltonian reads, 



dx 

o 7l' 



-ikx 



( x ) = <t>: 



k > 



HI 



^ fe | 2 + (l + fc 2 )|4| 2 



X] <^9 <V <^ 

9 9' 



k = — n k 



-q-q'-k 



Tik € Z 



(2.5) 



(2.6) 



The wavenumbers fc are dimensionless, the dimensionful momenta (conjugated to a; 1 ) are given by mk. 
The thermal average of any physical quantity = ©[</>, ir] in the canonical ensemble is written as 



(9[7T,0]) = 



JJ D(/>Dne- 0H ^e[4>,ir] 



(2.7) 



where J J DtfiDir stands for functional integration over the classical phase space and (3 = 1/T is the inverse (dimen- 
sionless) temperature in the dimensionless variables. In terms of the physical temperature, here defined as T p , the 
effective dimensionless temperature that will enter in the analysis that follows is given by 



T = — 



rn 



3 T P 



(2.8) 



As a consequence of the field redefinition available in the classical theory, the relevant variable for equilibrium ther- 
modynamics is T. Therefore for a fixed physical temperature T p we see from eq. ()2.8JI that the low temperature limit 
T <C 1 corresponds to the weak coupling limit and/or T p <C m. This will be relevant in the analysis below. 

Translation invariance (which is preserved by PBC) implies that averages of local observables O(x) which depend 
on 7r and <f> only at one point x, do not depend on x, that is (0(x)) = (©(0)) = (0). 

Furthermore, the fact that the Hamiltonian is the sum of a kinetic and a potential term, namely 

H[-K,4>]=r[K} + v[ ( j>] , 

entails that the average of observables represented by products of polynomials of the canonical momenta and the field, 
Q[4>, 7r] = 0i [4>] 02 M factorizes in the form 



with 



(6[<M> - (9iM>* (QaW), 



J Dcj> e-WM 



(e a W) w 



Moreover, since the it— integration in gaussian and ultralocal, it can be performed quite easily in most cases, leading 
to the purely configurational integral over <p = (f>(x) for < x < L. 



B. From quantum to classical field theory at equilibrium 



Before presenting our study of the non-equilibrium dynamics of the classical field theory, it is important to establish 
the equilibrium correspondence with the classical limit of quantum field theory. This will also allow us to make contact 
with non-perturbative treatments that are widely used in quantum field theory and that have a counterpart in the 
classical limit. We first consider free fields and then the non-perturbative Hartree approximation. 
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1. Free field theory 

The free-field expansion of the quantum field ip{x) and its canonical momentum Tl(x) for PBC in the interval 
< x 1 < L p is given by, 



tp(x) 
U(x) 



2uj„ 



a a e 



E 



a D e 



E fi ^°) ^ 



vV + - 



(2.9) 



(2.10) 



(2-11) 



The operators a q , aj obey canonical commutation relations. In free field theory in thermal equilibrium at temperature 



T p we have 



(a J a q ) = n q 



1 



?huJ q fTp ^ 



The classical limit is obtained for very large occupation number n q S> 1 which requires that ftw q <C T p . Therefore in 
the classical limit 



huj n 



(2.12) 



and we obtain the following thermal expectation values 



^-«) = 2^ [1 + 2 "" ] * f + rfi 



<n,n_,) 



hie, 



1 [1 + 2n 3 ] « T p 



(2.13) 
(2.14) 



In order to avoid the ultraviolet catastrophe in the classical statistical mechanics of the <fi 4 theory we introduce a 
ultraviolet momentum cutoff which we write as A p = irm/(2a), where 2a is a dimensionless lattice spacing. 
We then get the following results, 



(if 2 (x)) d ass = 

((p 4 (x)) claS s = 
(II 2 (x)) c i a ss = 



/ i' dip(x) \ « 

I 0x^ ' class 



IpST^-t f + °° 

Lp ^ w? - ^ 



dq 1 T p 
2-7T q 2 + m 2 2m ' 



3 Tl 

4 m 2 

If 
L 



p ^ lp J_> 2n 



dq _ T p 



\ ; -,l 2(1 p ' 



" t, 2^ ,,a ^ 1 v 



~ Ap dq q 2 
Y 2tt q 2 + m 2 2a p 2 



T p mTp 



(2.15) 



(2.16) 



where terms of the order 0{a) were neglected and we used Wick's theorem to compute (ip 4 (x)). The results (I2.15|) - 
(|2.16|l above lead to the classical result for the expectation value of the energy density in free field theory, 

1 T 

£ class = -j—(H mt \(n,(p)) dass = ~ ■ 

Furthermore, since the free theory is gaussian we find the ratio [ see eas. ^2.15l - ^2.16l ]. 



(y 2 (z)) 2 
(<p*(x)) 



1 



(2.17) 
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2. The Hartree approximation. 



The Hartree approximation is a very useful scheme that has been implemented as non-pert urbative treatment of 
the equilibrium as well as non-equilibrium aspects of nonlinear quantum field theories |l0ll2l| . This approximation is 
equivalent to a Gaussian ansatz for the wavefunctionals or the density matrix and leads to a self-consistent condition. 
Since it is one of the few non-perturbative techniques that are available and is a useful tool to study non-equilibrium 
phenomena, it is important to test its validity in the classical theory. Furthermore, we will discuss below exact results 
for the equilibrium situation in the 1 + 1 dimensional classical theory, which allows us to test the reliability of the 
Hartree approximation. 

In the Hartree approximation the quartic non-linearity is replaced by a quadratic form as follows 

^ 4 - *h V 2 { V 2 ) , (2.18) 

leading to the following effective mass and frequency (the momenta in the equations that follow are dimensionful) 

M 2 =m 2 + 3\((p 2 ) ; fl 2 q = q 2 + M 2 . (2.19) 

Thus the theory becomes Gaussian and the non-linearity emerges in terms of a self-consistency condition. The 
equilibrium occupation number and its classical limit is therefore given by 

and the Hartree self-consistency condition in the quantum theory is given by 

Using the classical limit en. (|2.2()|l and the result given by ea . (|2. 15fl we find the classical limit of the Hartree self- 
consistency condition 

(V 2 )H = ^- 1 (2.21) 

Assuming that in the high temperature limit (f 2 )H S> m 2 /A (to be confirmed self-consistently below) we find the 
solution of ea. H2.21|) to be given by 

^)H-T^(5) 2/3 = 0.436790...r^. (2.22) 
Since the Hartree approximation is a self-consistent Gaussian approximation, the ratio 

(4> 2 ? H i 



A4\„ 3 



H 



takes the same value than for free fields [ea. (|2.17(l ]. 
We also find the following results in the classical limit 



(H 2 (x)) H , cl = (2.23) 
2a p 



, f dipjx) \ . 

K \ dx 1 ) 



= {IL 2 {x)) H ,ci- M 2 (<p 2 ) H ,ci (2.24) 



with M 2 given in ea. l|2.19|l . In order to compare to the exact results obtained in sec. Ill Gl from the quantum transfer 
matrix, we summarize the results for the power spectra in the classical limit in the Hartree approximation in terms of 
the dimensionless field and canonical momentum [see eqs. (12.2(1 . I|2.6|l and (|2.8|l ] in the high temperature limit T 1, 

(\n q \ 2 )H,d = T 

T T 
<WVci = p + (|r) 2/3 = fc 2 + 1.310371... TV3 (2 - 25) 

These results will be compared to the exact results obtained from the transfer matrix in sec lll Gl and together with 
those provide a yardstick to establish criteria for thermalization and virialization. 
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C. Exact results from the quantum transfer matrix 



In this section we revert to the dimensionless variables given by eqs. (|2.2() - (|2.8|l . The expectation value of a 
7T— independent observable = 0[0] is given by the configurational integral 



(0(0)) = 



jD^e-OVM 0[ 



where 



V[<f>] = 



dx 



,-/3 v[ 



0' 2 + 2 + ^0 4 



(2.26) 



This is just an euclidean path integral for the path = 4>(x), < x < L, with PBC, so that the classical statistical 
average in eq. (|2.26l) can be interpreted as a quantum mechanical trace for the one dimensional anharmonic oscillator 
in imaginary time — ix with < x < L |22| . Namely, 



Tr (e- LH Xe[<t>]) 



-pv[, 



Tr e 



-LH 



(2.27) 



where is the family of (imaginary time) Heisenberg operators 



0(x) 



xH 



0(0) e 



-xH 



< x < L 



X stands for operator ordering along x. That is, 

X0(x)0(j/) = 6{x - y) 4>(x) 4>{y) + 6(y - x) 4>{y) j>{x) 
and TL is the quantum Hamiltonian 



nil)) i \ \ [m 



2 1 

+ 4 



0(0) 



Here 0(0) and f[(0) are canonical conjugate Schrodinger operators, that is [0(0), 11(0)] = i. 
It is convenient to make the following canonical change of variables, 



1 



Therefore [q,p] = i and 



<P(0) = VTq , Tl{Q) = —p 



(2.28) 



Hence, the calculation of classical statistical averages in eq. (|2.26(l reduces to quantum traces of the one-dimensional 
anharmonic oscillator in ea. H2.28|l with anharmonicity equal to the (dimensionless) temperature [see ea. l|2.8|l ]. 

The spectrum of H is an infinite set of discrete energy levels E n (T) associated with real eigenfunctions ip n (q;T) 
which solve the Schrodinger equation 



dq 2 



T 4 
2 q 



(q;T) = E n (T) ip n (q ; T) , n = 0,1,2,... 



and have defined parity: 



M-q-T) = (-i) n Mq;T) 



(2.29) 



(2.30) 



since TL is invariant under q — > —q, 

In the case of ultralocal observables 0[0(a;)], such as polynomials in <j){x), owing to traslational invariance (which 
translates into the cyclic property of quantum traces) we can express the thermal averages as single sums over diagonal 
matrix elements, 



(0(0)) 



1 



Z(T) 



-LE n {T) 



(2.31) 
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where 



with the proper normalization 



&mn{T) = I dq Q(VT q) ip m (q;T) ip n (q\T) 



+ OC 

dq [ip n (q;T)} 2 = 1 . 



For example, we get for 9 = [</>(a;)] 2 and 6 = [<^(a;)] 4 : 



(2.32) 



where, 

r+oo 



/ + CO 
dqq a iP m (q:T) ^ n (q:T) 
-oo 



The situation is more involved in the case of observables containing products of fields at different points. For instance 
the equal-time two-point correlation function takes the form 



' Tr[e- L *X4>(x)4>(y)] =-^^ (q nm f e -V>-\*-v\)E n <r)-\*-v\E m <X) (233) 



Z(T) 1 ^ '^ ya,i Z(T) 



Particular care is needed to treat local observables which contain powers of the field derivative 4>'(x), due to contact 
terms. Consider for instance [</>'(x)] 2 ; we have 

X [4>{x + dx) - 4>(x)} 2 = -i[4>(x), [H, 4>{x)] dx + ([H, 4>(x)} dx) 2 + ... = Tdx-T 2 U(x) 2 (dx) 2 + ... 

so that 

0'(x)] 2 =T{dx)- x -Te x ^p 2 e- X ^ 

The c— number singularity proportional to (dx)^ 1 reflects that the relevant configurations the ID functional integral 
are random walks with differentials proportional to (dx) -1 / 2 . Notice that the same singularity appears in the average 
of 7r 2 , since the tt— integration in gaussian and ultralocal, that is 

{n(x)n(y)) =T6(x-y) , (it 2 ) = T 5(0) = T (dx)^ 1 . (2.34) 

In terms of the lattice cutoff, 6(0) = 1/dx = l/(2a), therefore 

(O -L.IA. (2 .35) 



Since n = <j> 7 we obtain 



(<P 2 ) (4>' 2 ) = ^) Tr[e-^p 2 } = £ ^WCH > ( 2 ' 36 ) 



where, 



dq[i>' n {q;T)] 2 . (2.37) 

-OO 
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It follows from eas. (|2.36|l and l|2.34|) that 

(0 /2 > = f a -^f) E^ 2 ]«»( r ) e ' LEniT) • ( 2 - 38 ) 

For long rings (L ^> 1) the sums in eas. (|2.31|) - (|2.32l) arc dominated by the ground state considerably simplifying the 
results and leading to the following expressions 

(6) L ^ 9oo(T) + O ( e -WT)) ; l^i rp [q 2 ]oo(T) + ( e -L«-o(T)) 

(04) L^X t2 [qi]Qo{T) + ^-L^ ( ^ 2) _ ( ^ 2) L£l T [p 2 ]oo{T) + Q ( e -^o(T)^ 

where Wo(T) = E\{T) — Eq(T) is the first gap in the spectrum of Tl, namely the effective quantum Hamiltonian in 
the transfer matrix. Notice that, since all eigenvalues E n (T) as well as all differences co n (T) = £271+1 (T) — Eq(T) 
increase monotonically with T, these L ^> 1 approximations, as well as other to follow, are uniform in T. 



D. Classical and quantum virial theorems 

In equilibrium, the classical virial theorem is a consequence of the vanishing of the Boltzmann weight in the partition 
function for large coordinates and momenta. For the classical field theory under consideration, the classical virial 
theorem can be summarized as follows 

liWrn =T8 ij 6(x-y) (2.39) 
where i, j = 1,2 and 771 = <f>, 7/2 = 7T. More explicitly, 

(^0= T ^-^ = (^0 (^^))=°=(^^)>- (2 - 40) 

For example, 

V [ '5<P{y)/ JfDcbDne-^^ J J Dn e~P H ^<® { V) 

and integrating by parts in the last step. For the Hamiltonian (|2.4|) this implies, besides the already mentioned 
(Tr(x)n(y)) = T S(x — y) and the trivial ((j)(x)Tr(y)) — 0, the nonlinear relation 

(0(aO [-4>"(y) + <l>{y) + 3 ^)]) = TS(x - y) . 

Hence upon letting x —> y, taking the volume average and using the result (|2.34J) along with translation invariance, 
we find the classical thermal virial theorem 

(0 2 ) = {<j>' 2 ) + (0 2 ) + (0 4 ) . (2.41) 
When combined with the energy functional H[tt, </>] given in eq. 

J23JI, it yields 

M = f = (0-i(rt, (2.42, 

where E is the average energy in the canonical ensemble. Since N = L/(2a) is the total number of degrees of freedom 
on the lattice, we find that the temperature T is related to the (average) energy per degree of freedom as 

T=! + ^<0 4 >. (2-43) 

Therefore, if a (cf> 4 ) — > then the temperature is identified with the energy per site. We will study below the conditions 
under which the temperature can be directly identified with the energy per site. Comparing eq. 12.41fl to eq. I|2.36[) we 
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obtain 



T 



Tr [e- LH p 2 } 



As a consequence of eas. (|2.32[l . this entails the identity 

Tr[e- i7 ^ (p 2 -q 2 -Tq i )}=0 . 
This result is a direct consequence of the quantum virial theorem 

(p 2 ) nn = (qU'(q) )„„ 

for the case of our potential U(q) = ^q 2 + ^q 4 . 

While the virial identity eq. 12.41|) was obtained in thermal equilibrium, it is more general and it can be derived out 
of equilibrium by focusing on time averages. To see this, consider the following quantity: 



J(t) 



dx [4>(x,t)Y 



o 



Its second derivative with respect to time takes the form, 



lit) 



dx 



(x,t)(j)(x, t) + <j)(x,t) 



Which upon using the equation of motion (|2.5|) and after integrating by parts using the PBC yields, 



I(t) = / dx 
Jo 

The time average of a physical quantity is defined by 

e 



(2.44) 



lim — 

t^oo t 



dt' 0(i') 



For quantities which are total time derivatives of bounded functions this average obviously vanishes: 

— e(t) - e(o) 



de 

~dt 



lim 



t 







In the cf) 4 theory every classical trajectory is bounded for bounded initial conditions; therefore, taking the time average 
of ea. H2.44|) yields 



dx 







(2.45) 



This is the virial theorem for the </> 4 theory with PBC. If the time average is translational invariant in space, then 
the integral can be omitted and eq. I|2.45|l takes exactly the form of eq. (|2.41|l , with thermal averages replaced by time 
averages. Of course, the two averages should coincide if the ergodic hypothesis holds, but the validity of ea. (|2.45|) 
does not rely on this fact. This usually entails that the time scales for virialization are much smaller than those of 
thcrmalization [see sec. IIII Fl for a discussion on this for the present <f> A model]. 

Combining eq. (j2.45|l along with energy conservation yields the following alternative expression for the total energy 







E= [ dx 


Jo 





It should also be noticed that an infinite number of relations similar to ea. H2.45|l can be derived by considering the 
second derivative of quantities such as J Q dx [<f)(x,t)] n for n > 2. 

These results on virialization, which is different from the statement of thermalization, also provides a yardstick for 
the numerical evolution studied in section ITTll 
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E. Correlation functions and power spectra 

Let us now turn our attention to the equal-time correlation function l|2.33|) . which by translation invariance is a 
function only of x — y, that is (4>(x) (f>(y)) — G(x — y,T). Its Fourier transform 

G(k,T) = f dxe~ lkx G(x,T) 
Jo 

coincides with the thermal average of the so-called power spectrum of the field <fi, that is [see ea. H2.60 ] 

G(k,T) = (\4> k \ 2 ) . (2.46) 
Using completeness, J2 n ilmn) 2 = (q 2 )mm, we verify that 

i]TG(fc,T) = G(O,T) = (0 2 )(T). (2.47) 

k 

In the limit, L ^> 1 with x = O(L ), \L — x\ ^> 1, the correlation function G(x,T) can be approximated as [see 
eq. 



oc 

G(x, T)=TJ2 e - [B " (T) ^ Bo(T)I l x l (Mo„(T)) 2 + O ( e -[^m- B o(T)]L^ 

71=0 

= T £e-^W ([g] ,2n + i(T)) +o(e~^ L ) . 



(2.48) 



n=0 



where oj n (T) = E2 n +i(T) — Eq(T) and we used that (?o,2n = due to parity invariance [see ea. (|2.3U[> ]. Similarly, in 
Fourier space we find 

G^T) = £ ([ g ]o, 2 » +1 (T)) 2 + O (e^m) , (2.49) 

where k may now be treated as a continuous variable. 
Now the normalization eq. I|2.47|l is written as 

/ - G(k,T) = G(0,T) = TJ2 (Mo,2n+i(T)) = T[ g 2 ] o(T) = (0 2 ) . 

^ n=0 

This can be supplemented by two another sum rules related to ((f) 2 ) and to the (classical or quantum) virial identities, 
that is 

/ g e G{ki T ) = T £ u n {T) \— - u n {T)\ (M 0)2n+1 (T)) 2 = L - T [ P \ (T) = (0' 2 ) , 



2 1 2 

^a; n (T)([g] 0)2n+ i(T)) =- , £ c 2 (T) ([g] 0)2n+1 (T)) = [p 2 ] 00 (T) . (2.50) 

n=0 n=0 

Eas. H2.49|l - (|2.50(l imply that, for any temperature, G(k,T) vanishes for large k exactly as 



G(fc,T) fc =°°^, (2.51) 



just like the free field case. It should also be noticed that all the sum rules above are almost saturated by the first 
term in the sums, for any temperature. This appears evident from the data plotted in fig. ^ and will be shown in 
more detail in section IrTTH It follows that G(k,T) is very well approximated by the first term in its series eq. 
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The two-point correlation function of the conjugate momentum tt in the classical theory in equilibrium is given by 

(ir(x) tt(j/)) =TS(x - y) . 

which leads to a flat power spectrum for tt 

<I^-| 2 )=T, (2.52) 

This of course is a consequence of equipartition, and gives a criterion to identify the temperature: the height of the 
flat region in the power spectrum of tt. This identification will be very useful for the interpretation of the numerical 
analysis presented in sec. IIVI 



F. Low temperature limit 

In the low temperature limit T -C 1 the Hamiltonian ca. (|2.28(l becomes a harmonic oscillator with eigenvalues 
E n (0) = n + i and eigenfunctions 

_i 2 1 
ip n (q, 0) = c„ e 29 H n (q) , c n 



2«/2 7r V4 v ^! 



where the H n (q) are Hermite polynomials. As mentioned above, in the dimensionless variables the limit T <C 1 is 
equivalent to the weak coupling limit of the classical field theory in which the non-linearity can be neglected. 
Using the relation = VTq we find, for L ^> 1 and T — > 0, 

(0 2 ) = T [q 2 ] Q0 + 0(T 3 ) + O (e- L ) = | + 0(T 3 ) + O (e^) , 

(0 4 ) = T 2 [q% + 0(T 3 ) + O (e- L ) = ^ T 2 + G(T 3 ) + ( e - L ) , (2.53) 



and from eas. l(2~3l)|) - (|2~3T|) we find, 



2 ) = j + 0(T 3 ) + O (e- L ) 



Therefore, we recover the free-field theory result [see ea. (|2.17|l ] 

to^-i , h.T^5.'. (2.54) 

Similarly, for the matrix elements that define the two-point correlation function G(x, T) we find the following low 
temperature limit T — > 0, 

[ 9 ] 01 (0) = -L , [g] ,a, + i(0) = for Z>1, 
G(^T) = |e-H + 0(T 3 ) , G(fc, T) = p~~7 + 0(T 3 ) , (2.55) 



so that the free-field results 



are obtained. The result above for G(k, T) is the same as that obtained from the classical limit of the free field theory 
ea. fFWi after the rescaling fT^ - HT^i . 

We can evaluate the temperature as a function of the energy density E/L for a -C 1 using the low temperature 
formula ea. H2.53(l . We obtain from ea. (|2.43|l . 

T = 2af + f<^>. 
Using the results above, in the low temperature limit we find, 

It should be stressed that the low temperature limit is asymptotic, in the sense that the perturbation series in T has 
zero radius of convergence. So, also the low density limit E/L — > at fixed UV cutoff or the continuum limit at finite 
E/L are asymptotic expansions. 
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G. High temperature expansion 

For any T > the quartic terms dominates in the Hamiltonian 12.28fl and the quadratic term can be treated as a 
perturbation. In fact, we can make the following change of variables in ea. l|2.29[) 

/ 2 XV6 

so that ea. (|2.29|l becomes 



\ -^+e+(l) V \ 2 Xn^T-^)=e n (T-^) Xn ($:T-^) 1 n = 0,l,2,... (2.57) 

where 

E n {T) = r-} e „(r- 2 / 3 ) , ^„(? ; t) = (|J xn(C;T- 2 / 3 ) 

Since the new eigenvalues e n (T -2 / 3 ) and the new eigenfunctions Xn{£,]T~ 2 /' i ) are entire functions of T~ 2 / 3 , all 
equilibrium quantities multiplied by the appropriate power of T have convergent high temperature expansions in 
powers of T~ 2 / 3 . In particular, for large T we have to leading order [23) . 

/ j>\ i/ 3 f T\ 1 ^ 12 

E n (T)~l-\ e „(o), Mq,t)~ ( - J X n(y\0). 



d£ 2 ' " ' \T 



rr\ 1/3 /rpX 1/12 



Then the results of refs. |24j,|25j, such as 

e (0) = 0.530181045 . . . , ei(0) = 1.899836515 . . . , e 3 (0) = 3.727848969 .... 
prove to be very useful. We thus find for the relevant thermal averages, when T ^> 1, 

(0 2 ) ~ 2 1 / 3 T 2 / 3 [e 2 ]oo + O (V L TV3 ) (2.58) 

(0 4 ) * 2 2 / 3 T 4 / 3 K 4 ]oo + O (e- L Tl/3 ) , (0 2 ) - (0' 2 ) ^ 2- 1 / 3 T 4 / 3 [p 2 ] 00 + (e~ L T ' /3 ) . 

Solving numerically for the ground state xo (£ i 0) of the purely quartic oscillator and computing the appropiate integrals 
we find 

lim = 0.456119954..., 

lim ^rjz= lim X ' A =0.56104... 2.59) 
Thus we find following ratios for high temperature, 

(0 2 ) 2 (<i> 2 ) - (ft 2 ) , x 

lim ±f-L- = 0.37077 . . . , lim ^ ' lA )^ 1 = 1 . 2.60) 

to be the compared with their counterpart in the zero temperature limit [the second limit in eq. H2.60J) is a restatement 
of the virial theorem ea. (|2.41f) when ((f) 2 ) is subdominant] . 

The high temperature result for (tfi ) in ea. l|2.59|l allows to relate the energy density to the temperature using 
eq.(E32l as 



E ^ T 
1 ~ 2a 



1 - 0.28052- • - To 



3^ 3 



(2.61) 



Thus also in the high temperature limit T> 1 but with Ta 3 <C 1 the energy density is directly proportional to the 
temperature, E/L ~ T/(2a), signalling that the energy is dominated by the spacetime derivatives of the field. 
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For high temperatures, the correlation function G(x, T) takes the scaling form 

G(x,T)~2 1/3 T 2/3 g(xT 1/3 ) , G(fc,T) ~ 2 4/3 T 1/3 .9 (fcr~ 1/3 ) , 

where 

^)=£e-^W( yo , 2 m) 2 , g( P ) = ]>>^) 2 (2.62) 
with fi 2 /+i = 2~ 1/3 [e 2 i+i - e ] and 



2/o,2/+i = / y hi+i{y) <Po(y) dy . 

J — OO 

In particular, tlx = 1.087096267... and 0] y ,i = 0.600804942..., y ,3 = -0.032461289.... Since the matrix 
elements yo.21+1 decrease very fast for I > we can approximate eq. H2.62[> by the first term with the result, 

G(k,T) = (\fa\ 2 )~T 9 °- 988800 --- G^,^ ~T 2 / 3 0.4547894... e" 1 - 087096 - tV3 l-l . (2.63) 

v ' ' uv l ' k 2 + 1. 181779... T 2 / 3 v ; v ; 

Notice that in this approximation G(0,T) ~ 7 12 / 3 0.45478939 . . . while the correct high temperature limit ea. (|2.59|) 
is only 0.3% higher. Similarly, at large k the approximated G(fc,T) behaves like 0.9888 ... T/k 2 , while the exact 
expression is T/k 2 [see ea. (|2.51() ]. Taking into account the infinite radius of convergence of the perturbation series in 
T -2 / 3 and the fact that G(0, T = 0) reduces to the first term only, this first term approximation must be uniformly 
good in T. 

We can now establish the comparison between the exact (and approximate) results obtained above and those 
obtained in the Hartree approximation (in the classical limit sec lll B 2)l in the high temperature limit: 

(4> 2 ) ~ 0.45612... T 2/3 , (0 2 ) ff ~ 0.437... T 2/3 (2.64) 

^-- 0.37077... , 4-^- 0.333... 2.65) 

(</> ) (<P )h 

..7 |2 , 0.988800... T ~ 2 T 

{]M > "fc 2 + 1.181779... T 2 /3 ' (l ^ l)ff "fc 2 + 1.3103707...T 2 /3- {2M) 

Thus we see that there is agreement between the exact results and those in the Hartree approximation in the high 
temperature limit better than 5 — 10%. 

An important consequence from these results is that there is a strong renormalization of the frequencies. In the 
free field theory, which is also the low temperature limit, the frequency of oscillation is Uk = Vk 2 + 1, but in the high 
energy density, or high temperature limit, the effective frequency is uik = Vk 2 + 1.1818 . . . T 2 / 3 . This is an important 
observation because in a kinetic description the particle number or distribution has to be defined with respect to a 
definite frequency. It is noteworthy that the Hartree approximation does indeed capture very efficiently the frequency 
renormalization, thus suggesting that a particle number for a kinetic description can be defined with respect to the 
Hartree frequencies. 



H. Intermediate temperatures 



For generic temperatures only qualitative properties can be analytically established. For instance, all eigenvalues 
E n (T) as well as all differences w„(T) = £ , 2,i+i(T) — Eq(T) increase monotonically with T. This implies that the 
L>1 approximations of previous sections are uniform in T. 

To obtain quantitative expressions one has to resort to the numerical solution of the quantum anharmonic oscillator, 
which is nowadays an easy task on modern personal computers. Actually, accurate determinations of the eigenvalues 
for few values of the quartic coupling were obtained already in the late seventies (see for example refs.[3, EH)- For 
example, we have at T = 2 |24[ , 

E Q {2) = 0.69617582... , E x (2) = 2.324406352 .. . , E 2 {2) = 4.327524979 ... . 

In fig. ^ we plot in log-log scale our numerical determination of several equilibrium quantities. These plots nicely 
interpolate between the low temperature behaviour eas. H2.53fl - H2.54fl and the high temperature behaviour eas. H2.59fl - 
(f2.60fl . Notice that the ratio ((f> 2 ) 2 / ((j) 4 ) is monotonically increasing and only changes by 11% when T goes from zero 
to infinity [see eas. (|2.54fi . (|2.65fl ] , with more or less half of the variation concentrated at T < 2. 
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I. Expected dynamical evolution 

The exact results in the high temperature (large energy density) regime already suggest a preliminary picture for the 
dynamics in the case such that Ta 3 < 1. In order to present this preliminary picture, it is convenient to summarize 
the main exact results valid in the high temperature limit (and Ta 3 <C 1). 



T 
2a 

<- 2 > = f , <^> = f 
la la 



(tt 2 ) = — , ((j) 2 ) = 0.45612... T 2/3 , (0 4 ) = 0.56104... T 4/3 , ^f-f- = 0.37077 .. . (2.67) 

2a (<p 4 ) 

1 T 

1 - 1.122.. .(Ta 3 ) 1/3 J ~ — (2.68) 

i,~ ,o\ ,,-t 0.98880 . . . T . „ . 

<K| 2 > = T , (IM 2 ) = fc2 + L1818 ,.. r2/3 ( 2 - 69 ) 

where we have assumed Ta 3 <C 1 in ea. l|2.68|l . 

The important aspect gleaned from these exact results is that in equilibrium the nonlinear term is subdominant as 
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compared to the space and temporal derivatives. Namely, while the harmonic spatial and temporal derivative terms 
(tt 2 ) ~ (cj)' 2 ) ~ T/a ~ E/L, the interaction term (0 4 ) ~ T 4 / 3 , hence (^ 4 )/(tt 2 ) ~ (0 4 >/(^ 2 > - (Ta 3 ) 1 / 3 < 1. 

Consider solving the equation of motion (|2.5[1 with initial conditions (f>(x,t = 0), ir(x,t = 0) so that the Fourier 
transform of <j), tt have a power spectrum distributed on modes k <C A. Since only small wavevectors are excited 
in the initial state, this entails that for large energy density E/L 3> 1 it has to be that (f> 4 ~ E/L ~ T/a and 
a large fraction of the energy density is initially stored in the interaction term but a very small fraction is in the 
derivative term <p' 2 since the initial power spectrum for 4> 12 is localized at small wavevectors. Thus initially, <f> 4 ^> 4> 12 . 
Since the time evolution conserves the energy, the interaction between modes can only transfer the energy. If the 
time evolution leads to a thermal equilibrium state, then after or near the equilibration time scale it must be that 
4>' 2 ^> 4> 4 since in equilibrium the ensemble averages (which by the ergodic postulate are equivalent to time averages) 
are (</> 4 )/(0' 2 ) ~ (Ta 3 ) 1 / 3 < 1. 

The system will therefore evolve from the initial state to a final state of thermal equilibrium by transferring energy 
to modes of larger k, namely the initial power spectrum will necessarily broaden and energy will flow from the small 
k values to larger values in order to increase 4>' 2 . Interaction energy will be then transferred to gradient energy via 
a cascade of energy towards larger k, namely an ultraviolet cascade. As an equilibrium state is reached the power 
spectrum of tt, namely |7Tfc| 2 must become flat for all wavevectors up to the cutoff, since in equilibrium (|7Tfc| 2 ) = T 
[see ea. <|2.35[) ]. Since initially the power spectrum has been prepared to be localized at small values of the momenta 
the flattening of the spectrum must be a direct consequence of the ultraviolet cascade, namely more and more 
wavevectors are being excited by the interaction. Interaction energy is transferred to spatio-temporal gradient energy 
at the expense of the interaction term becoming smaller. We envisage this cascade towards the ultraviolet as a front 
in wavector space at a value k{t) that moves towards the cutoff as time evolves. For k <C k(t), | tt/c | 2 is approximately 
independent of k and eventually must saturate at the value T when the front k(t) reaches the cutoff. 

As argued above, initially (j> 4 3> <t> 12 since the initial power spectrum for the field is concentrated at wavevectors 
well below 1/a and the energy density (conserved in the evolution) is E/L ^> 1. However if the evolution leads to 
thermal equilibrium interaction energy is transferred to gradient energy via the ultraviolet cascade as described above, 
resulting in that <j> diminishes and 4>' 2 increases with time. Therefore there is a time scale to at which both terms 
are of the same order and the dynamics crosses over from being dominated by the interaction to being dominated by 
spatial derivatives. For t <^t , cj) 4 » <\> 12 while for t >• to, cp 4 <C 0' 2 . 

For t ^ to the theory is weakly coupled since the interaction term is much smaller than cf>' 2 and 7r 2 , the ensuing 
dynamics becomes slower, the non-equilibrium dynamics during stage is probably well described in terms of kinetic 
equations, since the small interaction guarantees a separation of time scales. 

These arguments that describe the expected dynamics are fairly robust and hinge upon very general features: i) 
energy conservation, ii) large energy density, iii) the initial condition determined by power spectra localized at small 
values of wavevectors, iv) the assumption that the dynamical evolution leads to a state of thermal equilibrium, v) the 
exact results obtained above in thermal equilibrium, which describe the final state of the dynamics. 

A detailed numerical analysis of the evolution presented below confirms the robust features of these arguments (and 
more) . 

J. Criteria for thermalization 

In the following section we undertake an exhaustive numerical study of the non-equilibrium evolution with the goal of 
understanding the dynamics that leads to thermalization and to assess the validity of the description presented above. 
We solve numerically the equation of motion l|2.5[) with given initial conditions. The initial conditions determine the 
energy density from which we can extract the equilibrium temperature either via ea. H2.56[) in the low temperature 
limit T < 1 or via eg. (|2. 61(1 in the high temperature limit T > 1. Since the classical field theory must always be 
understood with a fixed lattice (or UV) cutoff 2a, the equilibrium temperature is obtained by providing the initial 
conditions, which fix the (conserved) energy density E/L and the value of the cutoff a. Furthermore our primary 
interest is in studying the high temperature limit T > 1, hence large energy density E/L, and will always work 
with values of the temperature and cutoff such that Ta 3 <C 1 in which case eq. I|2.61|) entails that the energy density 
E/L ~ T/2a. 

In order to recognize when thermalization occurs we must define a consistent and stringent set of criteria. These 
are the following 

• The ensemble averages should give the same result as the temporal averages over a macroscopically long time 
plus spatial volume averages. That is, we check the validity of ergodicity. 

• The spatial (volume) and temporal average of ir 2 (x,t) should approach the canonical ensemble average (tt 2 ) = 
T/(2a) = E/L. 
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• The spatial (volume) and temporal average of 4>(x,t) , cj) 2 (x,t), [<f>'(x, i)] 2 , [cb(x,t)} 2 and [4>{x, i)] 4 should ap- 
proach their exact high temperature result eas. H2.59[) - i|2.59[l for high temperature and the behaviour plotted in 
figs. Hfor all temperatures. 

• The spatial (volume) and temporal average of the canonical momentum correlator must tend to T5(x — x'), 
which on the lattice would translate to T/(2a) [see eas. ()2.34[) - (|2.35[) ]. 

• The temporal average of |-7Tfc(t)| 2 (Fourier transform of the 7r-correlator) , which in equilibrium is given by 
(|^ fc | 2 ) = T. The temporal average of |0fe(i)| 2 (Fourier transform of the ^-correlator) must reach their thermal 
values. 



III. DYNAMICS OF THERM ALIZATION. 

Having provided an analysis of exact as well as approximate results of the cutoff classical field theory in equilibrium, 
we now pass on to the study of the dynamical evolution. In this section we study numerically the solution of the 
equation of motion (|2.5[1 with different initial conditions, of a given large energy density. By changing the initial 
conditions for a fixed energy density we are studying the evolution of different members of a microcanonical ensemble 
on a fixed energy (density) shell. 

Before embarking on the numerical study, we detail below our approach to solving the equations of motion by 
discretized dynamics in light-cone coordinates. 

A. Discretized dynamics in light— cone coordinates 

In order to solve numerically the evolution of the </> 4 theory it is necessary to discretize space and time. We choose 
to do that on a light-cone lattice 2(|. In this approach space and time are simultaneously discretized in light -cone 
coordinates with the same lattice spacing et-\/2~, thus preserving as much as possible of the original relativistic invariance 
of the field equation 

4> - <t>" + <f> + 4? = o. (3.1) 

Moreover, we choose a scheme where the discretized dynamics possesses an exactly conserved energy on the lattice. 
Given a space-time field configuration <f>{x, t), consider the two quantities 

1 2 1 2 

£±(x, t) =— \(f>(x + a, t) — (f>(x, t ± a)] + — \4>{x — a, t) — (f>(x, t ± a)] + 
2a 2a 

+ ^[l + (j) 2 (x,t±a)] [2 + <j) 2 (x + a.t) + <j) 2 (x - a,t)] -| . 

In the limit a — » 0, assuming <p to be smooth enough, we obtain 

£±=a (0 2 + |V0| 2 + 2 + i0 4 ) + 0{a 2 ) , 

identifying the leading term in a of both £ + and £_ with the energy density of the </> 4 theory. However, to higher 
orders in a they do differ; in fact 

£ + {x,t) -£_(M) = ^ + a )-^- a ) g( t) 

a 

where 

Q(x,t) = [(j>(x,t + a) + 4>(x,t- a)] |l+ ^[2 + 4> 2 (x + a,t) + 4> 2 {x - a,t)] J - <f>(x + a,t) - </>(x - a,t) . 
Hence, if Q(x,t) = 0, then £+(x, t) = £-(x,t) also on the lattice. In this case the total lattice energy 

E = ^£+{2na,t) (3.2) 

n 
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is exactly conserved in time, since it can also be written 

E = ^E-{2na + a,t + a) . (3.3) 

a 

This holds exactly on infinite space. If space is restricted to the segment [0, L] , suitable boundary conditions on 4>{x, t) 
are necessary; PBC are of this type if L = 2Na with N an integer. 

In conclusion, we may regard Q(x,t) — as a discrete field equation which conserves the total energy E. More 
explicitly, Q(x,t) = can be written as the recursion rule 

,/ , , n ,/ , , , (j>(x + a,t) + (j)(x - a,t) 

(p(x,t + a) = —<p(x,t - a) H — r (3.4) 

l + ^[2 + ^{x + a,t) + ^{x-a,t)} 

which evidently allows to propagate in time any configuration known in a time interval of width a. 

It is easy to check that Q(x, t) = indeed becomes ea. (|3.1f) in the continuum a — > limit. The order a is trivially 
satisfied, odd powers of a vanish identically as a consequence of the symmetry of ea. H3.4f) under a — * —a, while the 
order a 2 produces ea.|3.]|. 

Keeping up to 0(a 4 ) in ea. (|3.4|) yields, 

a 2 i~\ [<t>" ~ - (1 + 2 ) + ^{f" ~ ?) ~ 0W>' 2 + 4> + 0{a A ) . (3.5) 

To cast ea. (|3.4|l in a form suitable for numerical simulations, wc define the lattice fields F(n, s) and G(n, s) as 
F(n, s) = (f)(2na,2sa) , G(n, s) = 4>{[2n - 1] a, [2s + 1] a) , n,seZ. 
We then obtain the iterative system 

G(n,s) +G(n + l,s) 



F(n,s + 1) = -F(n,s) + 



1 + f [2 + G 2 (n,s)+G 2 (n + l,s)] 



„ x F(n-l,s + l)+F(n,s + l) ln 

G(n, s + 1) = -G(n, s) + 3 — ■ — . (3.6) 

V ' V ' 1 + f [2 + F 2 (n-l,s + l) + F 2 (n,s + l)] V ' 

with the PBC F(n + N,s) — F(n,s) and G(n + N,s) — G(n,s). As initial conditions we have to specify F(n,0) 
and G(n,0) for < n < N — 1. Once these values of the fields are specified, the iteration rules l|3.6|l uniquely define 
F(n, s) and G{n, s) for all s ^ 0. A comparison of this discretized dynamics with other more traditional numerical 
treatments of hyperbolic partial differential equations was performed in |27| . Here we only notice that this approach 
is particularly efficient, stable and accurate, specially when the continuum limit a — ► and very long evolution times 
are of interest. 

All observables of the continuum are rewritten on the lattice in terms of the basic fields F(n,s) and G(n,s). In 
particular, time and space derivatives are replaced by finite diferences. We always choose symmetric discretization 
rules such that lattice observables differ from their continuum limit by Ola 2 ). In the sequel, while referring to properly 
discretized observables, we shall keep using the continuum notation for simplicity, 



B. Initial Conditions 



We studied a variety of initial conditions with a fixed energy density in our calculations. In these studies the power 
is concentrated in the infrared; that is, 1 0fc(O)| 2 and | tt/c (0) | 2 are large for wavenumbers well below the cutoff A = j-. 
We considered the following sets of initial conditions: 

• Superpositions of plane waves: The initial fields have the form 

cj)(x, 0) = c; cos(fci x + 2 nji) , n(x, 0) = A dj cos(fej x + 2 7T 5j) . (3.7) 

i i 

The wavenumbers fcj = 2imi/L are chosen in the interval [0, fc max ] with fc max = 50tt/L — 25it/(N a) -C n/(2a). 
We will refer to these initial conditions as hard because these plane waves have sharp values for the wavenumbers. 
The power spectra \4>k\ 2 ', \^k\ 2 are sharply peaked at discrete (and small ) values of the momenta small compared 
with the cutoff A = 
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• Superpositions of localized wave packets: These are configurations of the form 

<j>(x, 0) = A^^ Ci w(x - Xi) , 7r(x, 0) = Ay dj w(x - yj) . (3.8) 

i i 

where we enforce the PBC by choosing, 

w { x ) = ^2 W o( k max[x -nL]) 
nGZ 

with wq(x) either the gaussian, wq(x) = e~ x , or the lorentzian, Wq(x) — . In practice, with our choice 
L r~j 20, only few terms in the sum are needed. These fields have support throughout Fourier space, but peaked 
as gaussians or simple exponentials at low wavenumbers k < fcmax ■ We will refer to these initial conditions as 
soft since these correspond to continuous and slowly varying power spectra. 

• Random on a fixed energy shell: A systematic way to decrease fluctuations is to average over initial 
conditions corresponding to a given value of the total energy. A microcanonical description gives equal a priori 
probability to all the configurations on the same energy shell. 

Thus, we choose smooth initial conditions as defined by eas. l|3.7p - H3.8fl and we average over them with unit weight 
choosing random values for the various coefficients. That is, we choose in ea. l|3.7|l the wavenumbers ki — 2nni/L 
at random both in number and in location, within the interval [0, fc max ] with fc max = 50n/L — 25ir/(Na) <C 
it /{2a). The positions Xi , in eci. H3.8|) are chosen at random in [0, L}. The phases in eo . (j^ . 7fl . ji ,Si, and the 
relative amplitudes Ci , di are also chosen at random in the interval [0, 1]. Finally, for any given realization of 
ki , 7i , Si , Ci , di, the overall amplitude A is fixed through the energy density E/L so all the configurations chosen 
are on the same energy shell. 

Typically, we performed averages over 30 initial conditions. 

For a fixed energy (density) the overall amplitude A in eas. l|3.7|l and l|3.8|> is fixed for any choice of Cj, di, ji, Si by 
the requirement that the initial configurations had always the same energy E. We considered several values of the 
energy density E/L ranging from 2.6 to ~ 10 4 . 

We notice also that all our initial configurations have vanishing total momentum J dx , which is a conserved 
quantity for PBC. 



C. Averaged observables 

The key observables in our investigation are the basic quantities 

0, 2 , 4 , 2 , , (3.9) 
as well as the power spectra of and tt, that is \(f>k(t)\ 2 and \nk(t)\ 2 , where, 

4>k{t) = j Q , = / ^e^7T(x,t), 

as in ea. l|2.6fl . 

In reference 1 13j the correlation function of the canonical momenta was studied which is particularly important since 
in equilibrium in the classical theory it is Gaussian. We here study a large number of independent observables and 
correlation functions since the criteria established above for thermalization based on the exact results apply to many 
different correlation functions. 

In accordance with the anticipation at the end of sec. IIII Al we are using here the continuum notation also for the 
Fourier transforms, although they actually are discrete Fourier transforms. 

The fluctuations of all these observables do not vanish upon time evolution. Hence for generic initial conditions 
they do not have a limit as t — > 00. These are fine-grained or microscopic observables. Typically there are several 
spatio-temporal scales, the microscopic scales correspond to very fast oscillations and short distance variations that 
are of no relevance to a thermodynamic description. We are interested in longer, macroscopic scales that describe the 
relaxation of observables towards a state of equilibrium. 

In particular the ergodic postulate states that ensemble averages must be identified with long time averages as well 
as spatial averages over macroscopic-sized regions. 
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To make contact between the time evolution and the thermal averages we need to properly average the microscopic 
fluctuations. 

First of all, for local quantities such as those in ea. (|3.9|l we take the spatial average. Secondly, we take suitable 
time averages of all key observables in the following way 

W(t) = - f* dt' - [ L dx[(f)(x,t')} 2 . (3.10) 

T Jt-T L Jo 

where r 3> a stands for the length of the time interval where we average. We find that the period of the fast 
(microscopic oscillations) themselves vary in time suggesting a sliding averaging in which r grows in time to compensate 
for the growth in the period of the fast time variation. We find that a practical and efficient manner to implement 
this averaging is to use 

T{t)=T(0)+Ct (3.11) 

where t(0); C small and positive with typical values ~ 0.1. In this way t 3> t(<) and the dependence on the initial 
values becomes negligible for practically accessible times. This method is quite effective in revealing general features 
of the (logarithmic) time evolution such as the presence of distinct stages characterized by well separated macroscopic 
time scales. 




Altogether, we denote the results of all coarse-grainings simply with an overbar (not to be confused with the 
complete time average of section III D|l to avoid cluttering of notation. For example we have in the case of random 
initial conditions on the energy shell 

where the superscript W labels the M different choices of smooth initial conditions, all with the same energy density 
E/L, we have taken M ~ 30 in our study. Likewise, 

1 *f r t 

N 2 (*)= mE ~ / dt '\^\t')\ 2 (3-12) 
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log t 

FIG. 3: 4?{t) as a function of the logarithm of time for T = 69.2, E/L = 3424, L = 20 and a = 0.01. 



where we used the reality condition 7Tfc = Trl k . Analogous expressions hold for <f>, <fi 2 , (f> 4 , ir 2 , 4>' 2 and J 



In particular, due to the linearity of these averages, we have the sum rules: 



+ A 



dk -~ 



\M 2 (t) = 02( t ) 



+A 



dk 



-A 2n 

Where UV cutoff in the light-cone lattice is A = n/(2a) 



dk 



\t k \ 2 (t) = n 2 (t) , / ^k 2 \4> k \ 2 (t) = <b' 2 (t). (3.13) 

A Zn J-A Zn 



The power spectra |0fe| 2 (i) and |7Tfc| 2 (£) are connected to equal-time correlation functions of <f>(x,t) and n(x,t), 
respectively, just as it happens at thermal equilibrium [see e.g. ea. l[2.46)l ]. We have, for instance 



\<t> k \ 2 {t) = dxe 
Jo 

M 



-ikx 



X,t) 



X. 



(3.14) 



with similar relations between |7Tfe| 2 (i) and or irir(x,t) and between k 2 \<fi k \ 2 (t) and (p'<p'(x,t). 
It is useful to define the normalized power spectrum of it, 



P(k,t) 



\^\ 2 (t) 

n 2 (t) 



which describes the distribution of power over the wavenumbers, it is normalized so that, 

r+A 



f dk „ , 
/ — P(k,t =1 
J-A 2n 



(3.15) 



(3.16) 



Another important quantity of paramount importance to describe the cascade described in section Til II above is the 
average wavenumber 



k(t) 



+A dk 

£ \k\P(k,t) 

-A IK 



(3.17) 
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In thermal equilibrium (|7Tfe| 2 ) = T, (ir 2 ) — T/(2a) — TA/tt, therefore in equilibrium k eq = A/2. 

The physical significance of this quantity becomes obvious by considering the situation in which the power spectrum 
\nk\ 2 is approximately flat in a region of wavevectors |fc| < ku and negligible elsewhere, namely |7ffc| 2 cx Q{k 2 M — k 2 ). 
In this case k(t) = The relevance of this effective quantity will become clear below when we study in detail the 
process of cascade of energy towards the UV, described in section lTTTl where it will become clear that 2k(t) determines 
the front of the ultraviolet cascade as described in section lrTTI 

To summarize, we perform spatio-temporal averages in the isolated system with fixed energy (density) which is 
equivalent to microcanonical ensemble averages at long time by the ergodic postulate. In the thermodynamic limit, it is 
expected that microcanonical and canonical ensembles will yield the same equilibrium results provided the equilibrium 
temperature T is identified with the energy density as per eas. (|2.56[l or l|2.61|l in the low or high temperature limit 
respectively. 

While in refs. [l3| and |27| only the variance of tt 2 was studied as a measure of gaussianity and thermalization, 
we study many different correlators, since as described above the exact solution in thermal equilibrium furnishes a 
stringent set of criteria for thermalization. 

D. Time evolution of basic observables 

We use the lattice field equations, eq. l|3.6[> . to evolve the initial configurations in time and compute the time and 
space average of the basic quantities (|3.9|l as a function of time for 10 < t < 10 s and 0.1 < a < 0.0001. 
We confirmed that the lattice energy eas. (|3.2(l - (|3.3|l is indeed conserved with large accuracy. 

Figs. El El display 4>' 2 (t), ^(t), and ^(t), respectively, as functions of time for E/L = 3424, L = 20 and 

a = 0.01 corresponding to an equilibrium temperature T ~ 68.5. Figs. [3] and [51 display the same quantities as 
functions of time for E/L = 50.28, L = 20 and a = 0.01 

In figs. I2l5l the initial conditions are the plane waves eq. l|3.7[l with parameters that fix tt 2 ~ </> 4 ~ E/L initially. No 
average over initial conditions is performed and we used the sliding time average with linearly growing time intervals 
as in ea- H3.11|) . 

These figures clearly reveal the expected dynamics as described in section lTTTI above. Initially <fi' 2 (t) is small reflecting 
the fact that the initial conditions determine a power spectrum localized at wavectors k <C \/a. The mode mixing 
entailed by the interaction is transferring power to larger wavevectors, thus effectively transferring energy from the 
interaction term, which diminishes, to the spatial gradient term which increases. As is clear from these figures, all 
magnitudes tend to a limit for late times. The late time limits are the thermal equilibrium values, as we discuss below 
in detail, 

The growth of 4>' 2 {f) at the expense of the interaction term </> 4 (i) shows that thermalization is a result of the flow 
of energy towards larger k modes, namely the ultraviolet cascade ultimately leads to the thermal equilibrium state. 

We find three distinct stages of evolution for a wide choice of initial conditions. 

• A first stage with relatively important fluctuations and whose precise structure depends on the initial conditions. 
Such stage can be seen in figs. 12151 for lnt < 5, that is t < 200. The value ~ 200 turns to be independent of 
both a and E/L as long E/L is not very small. Namely, we find U ~ 200 for 0.1 > a > 0.0002 and E/L > 10. 
For lower values of E/L, U increases sharply as shown in figs. I""l and 12*51 In addition, for such low values 
of E/L, ti becomes dependent on the details of the initial conditions but is not relevant for our study which 
focused on the high density case. During this first stage the transfer of energies via the cascade begins to be 
operative and becomes most effective at the time scale ti. We see that U can be identified with the time scale 
at which the interaction and gradient terms are of the same order and the crossover from a strongly interacting 
<f> 4 ^> <\>' 2 to a weakly interacting theory occurs. Thus the first stage corresponds to t < ti ~ 200 during which 
the cascade is the result of large interaction energy which is redistributed to larger wavevectors via mode mixing 
and the dynamics is dominated by the interaction. 

• There is a second stage where (j)' 2 (t), 4> 2 {f), 4> 2 {f) and 4> 4 {t) are about the same order ~ E/L. During this 
stage there is a crossover from a strongly to a weakly interacting theory since at the end of this second stage the 
spatio-temporal gradient terms are much larger than the non-linear term. During this second stage the cascade 
is very efficient in redistributing the energy. 

In figs. 12151 this second stage corresponds to 5 < lnt < 10 for a = 0.01. 

We find that the behaviour of the observables during this second transient stage depends to some extent on the 
initial conditions. For hard initial conditions [as ea. l|3.7|l ] we find much steeper curves for 4>' 2 (t), (f> 2 {t), 4> 2 {f) 
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and 4> 4 (t) than for soft initial conditions [as eas. (|3.8(l ]. Thus this second stage corresponds to an interval after 
the time at which the spatial gradient and interaction terms cross. 

While the details of the dynamics during this stage depend on the initial conditions, the presence of this stage 
during a time interval 5 < lnt < 10 ~ Into is fairly robust. We found such intermediate stage for all types of 
initial conditions, and a wide range of lattice spacings and energy density. 

During this second stage 7r 2 (i) varies with time more slowly than <p' 2 (t), 4> 2 (t) and 4> A {t). 

• After this transient second stage there is a third stage where these physical quantities approach their thermal 
equilibrium values. The third stage starts at a time scale to ~ 50000. The value of to turns to be independent 
of both a and E/L as long E/L is not very small. That is, for E/L > 10. 

In particular, Tt 2 (t) slowly decreases to its asymptotic value which agrees with the thermal equilibrium value 
T/(2a), that is, using ea. l(235|l : 



lim ir 2 (t) 



T 
2a 



E 
L 



+ 0{a 2 ) 



(3.18) 



The dynamics during this third stage is also driven by the ultraviolet cascade but unlike the earlier two stages 
wherein the details depend on the initial conditions, we find that this third stage is described by a universal 
cascade independent of the initial conditions as we shall see in the next section. In particular, the results 
obtained during the third stage for sharp initial conditions [ea. (|3.7ll ] are independent of the chosen wavenumbres 
ki provided they do not approach the lattice cutoff. The third stage ends at a time ~ t\ when the lattice size 
effects start to play a role. For the cases depicted in figs. 12131 and figs. 14151 we have lnii ~ 12 and lntj ~ 17, 
respectively. 

• The fourth and final stage extends for times beyond t\. Thermalization is reached here, strictly speaking, only 
for infinite times. Since the dynamics feels here the details of the discretization adopted, we shall not study this 
stage in detail. 
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FIG. 6: </> 2 (t) as a function of the logarithm of the time t for E/L = 3.67, L — 20 and a = 0.01. For such a low value of E/L 
the cascade starts much later for logio ~ 17 



Figs. 17181 display (f>(t) and log \<fi(t)\ as a function of the logarithm of the time for E/L — 3424 and E/L — 50.28 for 
L = 20 and a = 0.01, respectively. 
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We see that the relaxation of <f> towards its thermal equilibrium value (4> — 0) is different from the other physical 
quantities previously discussed. This is due to the fact that the vanishing of <f> is connected to a symmetry of the 
model. We find that <j>(t) relaxes as l/t for t < t (during the first two stages of thermalization) and as ~ \j\ft for 
t > to (during the two subsequent stages). 
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FIG. 7: 0(t) and log |0(t)| as a function of the logarithm of the time t for E/L = 3424, L = 20 and a = 0.01. 



E. Time Evolution of the Correlation Functions 



We now turn our attention to the study of correlation functions. 

According to the Fourier transform relationship ea. H3.14|) between the power spectrum |0| 2 (fc, t) and the equal-time 
correlation function (jxf>(x,t), there are two approaches to the numerical evolution of such quantities (we specialize 
here on <f> but the discussions applies equally well to tt). We extract the field <fi(x,t) from the lattice fields F(n,s) 
and G(n, s), Fourier-transform it to <pk(t) and then perform all needed averages on |</>fc(i)| 2 . Or we directly compute 
averages of the correlations of F(n,s) and G(n,s) and extract from them the correlations <fc<j>(x,t) and WW(x,t). We 
found that both methods yield the same numerical results (see Appendix). 

Moreover, when using the approach with growing time averages as in ea. (|3. 1 ljl with a unique initial condition, one 
realizes that the simply time-averaged correlation 

1 l* 

- / dt' (j>(x,if) 4>(x',t') . (3.19) 

T Jt-r 

very soon (in the logarithm of time) becomes translation invariant, that is a function only on the distance \x — x'\, 
making the time consuming space average unnecessary. 

We plot in figs. l9l HtH and ITT1 thee example of the typical profiles of (fxj>(x,t), fc 2 \<pk\ 2 {t) and \iTk\ 2 (t), respectively, 
for one given choice of parameters. We see from figHJ] that the angle point at x = 0, characteristic of the equilibrium 
correlation function G(x) [see eas. (|2.33|l and ()2.48(l ] is developping. This appears more evident in Fourier space, fig 

I1UI since the region where k 2 \<fik\ 2 (t) is large and almost constant is spreading towards the UV cutoff. Likewise we 
see the same UV cascade for |7Tfc| 2 (£) in fig^] The power spectrum of the canonical momentum w will be discussed 
in detail in the next section. 
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FIG. 8: <j)(t) and log \4>(t)\ as a function of the logarithm of the time t for E/L — 50.28, L — 20 and a = 0.01. 



F. Early Virialization 



We depict in fig. El the quantity, 



This quantity vanishes when the virial theorem is fulfilled [see ea. (|2.41fl ]. It turns out to be negative for finite times 
and nonzero a. We see from fig. ^Jthat |A(t)| starts to decrease at times earlier than to. Therefore, the model starts 
to virialize before it starts to thermalize. |A(t)| keeps decreasing with time and tends to a nonzero value which is 
of the order C(a 3 ) for t — > oo. This is to be expected since ea. H2.41|l only holds in the continuum limit and receives 
corrections in the lattice. 

IV. THE ENERGY CASCADE 

We describe here the flow of energy towards higher frequencies leading towards thermalization. Such cascade turns 
to be universal (independent of the lattice spacing and of E/L) and exhibits scaling properties within a wide range 
of time. 



A. Power spectrum of n and the universal cascade 



The chosen initial conditions eas. (|3.7fl - H3.8|) are such that the power is concentrated in long wavelength modes k 

well below the ultraviolet cutoff A = it /(2a). Therefore, |7Tfc| 2 (0) ; |</>fc| 2 (0) are concentrated on small k. During the 
time evolution the non-linearity gradually transfers energy off to higher fc— modes leading to the ultraviolet cascade 
as discussed above. 

The typical behaviour of | vr/c | 2 (i ) is shown in figs. 1131 1141 and 1151 In these figures we have averaged over time with 
constant intervals r, and on the initial conditions as defined by ea. 13. 12(1 . This plots have also been smoothed by a 
moving average over wavenumbers with averaging intervals of size 327r/L [see the Appendix for details]. 

The time evolution of the power spectrum |7r ( t| 2 (t) features the ultraviolet cascade in a very clear manner. The 
gradual transfer of energy to larger wavevectors results in a the formation of a central plateau which spreads over 
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FIG. 9: The correlation <fxj>(x,t) as a function of x at four different times. It takes the thermal form eg. 1)2. 630 for an effective 
temperature T e s(t). 



higher wavenumbers as time grows decreasing its height. This plateau ends abruptly at a value of the wavector which 
determines the front of the cascade. From the definition of k(t) given by ea. 1)3. 17)1 and the discussion that follows it, 
it is clear that this cascade front is given by 2 k(t). 

The limiting form of | tta; | 2 C f° r t ~* 00 i s as expected from thermal equilibrium, ea. 1)2.52)1 [see fig. 115) : 



lim 



T 



(4.1) 



The second and third stages described in sec. III.D for t to correspond to the steady flow of energy towards 
higher wavenumbers with a steady increase of the average wavenumber k(t) towards its asymptotic value of thermal 
equilibrium [see figs. E] an d decrease of the height of the plateau. Thus the front of the cascade 2 k(t) advances 
towards its limiting value, given by the cutoff 2 fc(oo) = A = j- leaving behind a wake in local thermal equilibrium 
at a temperature T e g(t) corresponding to the height of the plateau. 

Figs. 1131151 clearly show that there are only two relevant wavevector scales in the power spectrum |7Tfc| 2 (i): the 
front of the cascade 2 k(t) and the momentum cutoff A. Neglecting fluctuations in the plateau, it is clear from these 
figures that for k < 2 k(t) the power spectrum is flat just as in the equilibrium case, and for 2 k(t) <C A the shape of 
the power spectrum is very simple and insensitive to the cutoff. 

We therefore reexpress the distribution P(k,t) defined by eg. ()3. 15)1 in terms of dimensionless ratios as follows, 



P(k,t) 



k(t) 



F(u, v) , u 



k 
W) 



v(t) 



2k(t) 
A 



4 - 

- kit) a 

7T 



(4.2) 



with the asymptotic behavior of v(t) given by 



v(t) 1 =°° 1 

Our numerical analysis shows that during a fairly large time interval t\ 3> t ^> io ; during which 2 k(t) <C A, the shape 
function F(u, v) becomes universal, namely independent on the initial conditions, including the energy density E/L 
[see fig. an( l independent on the UV cutoff as well [see fig El- The new time scale ti is determined by when the 
front of the cascade begins to reach the cutoff, namely 2 fc(^i) ~ A. At this point the wake in the power spectrum 
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behind the front of the cascade corresponds to the plateau at the temperature T, namely T c g(ti) ~ T. The time scale 
t\ marks the end of the third stage. Beyond it, nonuniversal (a-dependent) effects become relevant. 

Our exhaustive numerical analysis leads us to conclude the following picture for the cascade and the process of 
thermalization. For a wide range of initial conditions, lattice cutoff and energy density, there exists a scaling window 
for the cascade, namely a time interval t\ 3> t 3> to characterized by the following properties: 



For ti ^> t to so that 2 k{t) C Awe find that k{t) is described by power-like behaviour in time, as shown for 
instance in fig lKil namely 

k(t) ~hit a ■ (4.3) 

where a is slowly growing with E/L, typically 0.21 < a < 0.25 for E/L > 10 , while hi ~ (E/L) 7 , with 
7 ~ 0.25, for E/L > 10. The numerical evidence suggests that both a and 7 depend on the initial conditions 
only through E/L. 

• During the interval t\ ^> t ^> to the front of the cascade is far away from the cutoff, namely u«l and F(u, v) 

depends very weakly on v therefore F(u, 2^^-) is almost time independent. Hence, during this interval we 
can make the approximation F(u, v) ~ F(u, 0) so that P(k, t) satisfies the following scaling law with great 
accuracy [see fig. El 

Hence P(k,t) satisfies a renormalization group-like relation for its time dependence, namely 



and F(u, 0) plays the role of scaling function. 
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FIG. 11: | -TTfe | 2 ) vs. k at the same times of figs. ISland Hni It is described in the central plateau by the flat thermal spectrum 
ea, 12.52| lfor an effective temperature T e a(i). The height of the plateau is identical to fie llOl as expected. 



The behavior of k(t) as a function of the logarithm of time for several values of the lattice cutoff a displayed in 
fig llOl shows a power-law behavior given by ea. (|4.3|l in the window t% t 3> to and a saturation for times larger than 
t\ (which depends on the lattice cutoff). 

Proposing that, for very small lattice spacings and late times, a k(t) is solely a function of the combination at a , 
namely 



a k(t) = h(s) , s = a t a , a 
h(s) ~ hi s , s — > ; 



, t 



(4.6) 



we find numerically the function h(s) displayed in fig. EDfor E/L = 5200 and in fig. |2]]for E/L = 520 and E/L = 104. 
These results show quite clearly that a does depend, although quite weakly, on E/L. Actually, as explained in the 
Appendix, this scaling-based approach is more effective in the determination of a than direct fitting. Moreover, these 
results show that h(s) is not universal, but depends on the initial conditions at least (and most likely only) through 
E/L. 

In fact, it turns out that the data for h(s) all collapse very well on a unique profile upon rescaling s — > [E/L) 1 s, 
if 7 = 0.25. In other words, for the average wavenumber the following double scaling form holds true when E/L is 
large enough 



k(t; a; E/L) ~ a 



-'*(•'■ (!)' 



(4.7) 



where h(s) is now a universal function of order 1 with h(s) ~ h\s as s — + and hi ~ 1.1 [our numerical reconstruction of 
h(s) is reported in fig. |^. Therefore, a fairly good approximation for k(t) during the universal cascade (ii ^> t 3> to) 
is given by, 



k(t) 



1.1 



(4- 



with 0.21 < a < 0.25 according to the value of E/L [see Table I]. 
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FIG. 12: Virialization: the logarithm of x |A(t)| vs. the logarithm of the time for E/L = 50.3, 728 and 3424 with 

a = 0.01 and for E/L = 34861 with a = 0.001.° 

The function h(s) saturates for infinite time to the value it /4 which translates into the maximum value for the front 
wavevector k max = ir/Aa — A/2. Actually h(s) in practice saturates already when > tt/2 (see fig. I21|l . Therefore we 
extract the thermalization time scale 



ii ~ A 1/a (E/L)-' ,ia , 



(4.9) 



at which 2fc(ii) < A. For t ^> ii the power spectrum |7ffc| 2 features a plateau for all wavevectors up to the cutoff, 
thus describing the thermal equilibrium state. 

As far as the shape function F(u, v) is concerned, we see from fig. ^|that for v <C 1 it features a bell-shaped profile 
with exponentially small tails for large u, 

F(u, v)~C exp[-7(«)|«|] , |u| > 1 , 

with 7(0) ~ 2 and C a constant. As v increases, the tails decrease even faster, the lateral walls steepen and the central 
region flattens. Finally, as v — ► 1, which corresponds to complete thermalization, we have [see eqs. (0~TJ) and (jTVl} ] 



lim P(k, t) = 2a , F(u,l) 

t—>oo 



tt/2 \u\ < 2 
\u\>2 



(4.10) 



to be compared to the approximation F(u, v) ~ F(u, 0) valid in the scaling window. Strictly speaking, the function 
F(u,0) is truly universal in the continuous limit where v = — k(t) a — > 0. However, one can see a very good 
approximation to F(u, 0) in fig. ^|for the plot at the smallest a — 0.000625. 

Finally, by construction, F(u,v) is even in u and must satisfy, following eas. l|3.16|l and (|3.17|) . 

2/v p2/v 

duF(u,v)= / duu F(u,v) = n , 
Jo 

for any value of v. 

Summarizing, after the transient stage, the effects of the initial conditions on P(k,t) are entirely accounted for 
by the average wavenumber k(t), which fixes the scale of the universal ultraviolet cascade with shape described by 
F(u, v). For v <C 1, namely when the front of the cascade is far away from the cutoff, the function F(u, 0) is universal 
as shown explicitly by the collapse of the data for several values of the initial energy density onto one function in fig. 
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FIG. 13: |^k| 2 (i) cascade still far from the UV cutoff. 



B. Time dependent effective temperature and local equilibrium 



800 



The shape and universality of the power spectrum |7Tfc| 2 (t) featuring a flat plateau behind the wake of the cascade 
leads to the definition of the effective temperature T e g (t) as the height of the plateau. The interpretation of this 
temperature is that for t 3> to but well before the true thermalization time t\ , there is an description in terms of local 
equilibrium (local in Fourier space) at the effective temperature T c g(t). 

In order to firmly state this conclusion, however, we must understand if the other criteria for thermalization 
established in section ITTjl are fulfilled. 

Combining eas. (|3.15j) and l|4.2f> we can write the tt power spectrum as 



1 1 ^ ' ' k(t) Vjb(t)' A J 
For to <C t <C ti, within the scaling window, this is well approximated using ea. (|4.3ll as 

A: 



ifpta t ) ~ ze!W f (— 







Now, for any fixed k in the bulk of the cascade and large enough t (but still t -C ii), taking into account the flat 
profile of F(u, 0) around u — 0, we can write 



hx t a 



F(0, 0) = Teff(t) , 



which is indeed k— independent. Moreover, the observed change in time of the integrated power spectrum 7r 2 (i) is not 
as important [see figs. I^Jand^]. It slowly decreases to its asymptotic value which agrees with the thermal equilibrium 
value T/(2a) [see ea. H3.18j l]. Hence, to leading order in t (but always with t <C ii) we may write: 



Teff(t) 



E F(0,0) 
T hi t a 



(4.11) 
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FIG. 14: |7ffc| 2 (t) cascade reaching the UV cutoff. 



Finally, using the large E/L behavior of hi [see ea. H4.7f) ]. we arrive at 

Teff^Trd) 1 ^-", (4.12) 

The constant T\ = F(0,Q)/hi is estimated to be T% ~ 1.55, very close to ir/2. We recall also that a is weakly 
dependent on E/L for E/L > 10, saturating approximately at 0.25 for large E/L, and that 1 — 7 ~ 3/4 quite 
precisely from the scaling argument of ea. H4.7fl . Therefore, fairly good approximations for t\ and T c g(t) are given by, 



The identification of the effective temperature allows us to study whether the other observables, such as 2 (i), 
4> 4 (t) or the correlation function <jxj)(x,t), depend on this effective temperature in a manner consistent with local 
equilibrium. A numerical analysis reveals that indeed to high accuracy we have [see figs. 1221 and 12*5 ] 

^(t) = <^>(T eff (i)) , ^(t) = (0 4 )(T cff (i)) , (4.14) 

in terms of the equilibrium functions (0 2 )(T) and (^ 4 )(T) plotted in fig. early times (but t ^> t ), that is high 
effective temperature these reproduce the thermal equilibrium results given by eq. (f2.67f) , while for late times and large 
UV cutoffs, that is low effective temperature, these agree with the free field results (f2.53fl - (f2.54fl . 

A more precise analysis can be performed on (fxf>(x,t), or better on its Fourier transform \<p\ 2 (k,t). Indeed, as 
evident from fig. |3| we find that for k <C 2 k(t) -C A and t > t 

\$(k,t)~G(k,T eff (t)), (4.15) 

in terms of the equilibrium distribution G(k,T). For k too close to the front of the cascade k < 2k(t) effective 
thcrmalization has not occurred yet, but this does not affect significantly quantities like 4> 2 (t) in which all wavenumbers 
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FIG. 15: |7t"fc| 2 (t) cascade close to complete thermalization. 

are summed over, since the equilibrium G(k,T) vanishes as T/k 2 for large k. It instead affects significantly quantities 
like <f>' 2 (t), since the <j> equilibrium power spectrum k 2 G(k, T) goes to T for large k. 

The physical meaning and interpretation of the effective temperature and the description in terms of local equilib- 
rium at this temperature behind the front of the ultraviolet cascade applies solely to long wavelength physics. Namely 
the wavevectors behind the front of the cascade k < 2 k(t) can be considered to be thermalized by the mode mixing 
entailed by the interaction. Physical observables that do not depend on the cutoff are described in terms of this local 
thermal equilibrium concept. 

This description in terms of a universal cascade in local thermal equilibrium at an effective temperature T c g(t) 
applies for t ^> to up to the time t\ at which the front of the cascade 2 k(t) reaches the cutoff. From this time onwards 
(during the fourth stage, see sec. IIII Dfl the power spectrum (|7r/c| 2 ) is no longer universal and is sensitive to the cutoff. 

There is a clear separation between the time scale to ~ 50000 where the cascade forms with the effective time- 
dependent temperature T e g (t) and the much longer time scale t\ oc a~ 4 which signals that the front of the cascade is 
near the cutoff and the end of the universal cascade. 

Thermalization does continue beyond t\ during the fourth stage [see sec. IIII D| and the power spectrum (|7Tfc| 2 ) 
eventually becomes T 0(A — k) with T the true equilibrium value of the temperature (determined by the energy 
density) at infinite time. Thus true thermalization takes an infinitely long time. 

Hence the dynamics for long-wavelength phenomena can be described to be in local thermodynamic equilibrium 
for t > to at a time dependent temperature T e g (t) while short wavelength phenomena on the scale of the cutoff will 
reach thermalization only at a much later time scale t^>t\. 



limit 0.37077 . . . [eq 



The ratio plotted in thermal equilibrium as a function of T in fig. ^ provides a simple test of effective 

thermalization. This ratio increases monotonically from its zero temperature value | up to its infinite temperature 

Therefore, a necessary condition for effective thermalization is that 



1 < < 0.37077. 

3" 04( f) " 



(4.16) 



We depict in fig|2S]the ratio %A (t) as a function of the logarithm of the time t for three values of E/L. Effective 
thermalization may start only when the ratio falls within the inequality l|4.16|l . We see from fig. [5S]that for E/L > 10 
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FIG. 16: Numerical evaluation of the average wavenumber k(t) for different values of the UV cutoff, at fixed energy density. 
When a = 0.01 the system almost thermalizes [see fig. 115) and k(t) saturates to its equilibrium value A/2. When a = 0.005 
the bulk of the cascade reaches the UV cutoff, causing a significant bending. For smaller a the bulk of the cascade is still far 
away from the cutoff at the times depicted here. 



this happens around hit ~ 9 which is about the begining of the universal cascade [see sec. HTTP] . For E/L = 3.669 
effective thermalization starts much later around lni ~ 12. More generally, we find for E/L < 10 that thermalization 
is significantly delayed. 



V. CONCLUSIONS AND DISCUSSIONS 



In this article we have studied the approach to equilibrium in the classical tfi 4 theory in 1 + 1 dimensions with the 
goal to understand the physical process that lead to thermalization and to establish criteria for the identification of 
a thermalized state in a strongly interacting theory. After discussing the classical theory as the limit of the quantum 
field theory, we exploited the equivalence of the classical partition function to the transfer matrix for the quantum 
anharmonic oscillator. A body of established results on the spectrum of the quantum anharmonic oscillator allows us 
to obtain exact results for low, high and intermediate temperatures which furnish a yardstick and a set of criteria to 
recognize thermalization from the physical observables. We compared these exact results to those obtained for the 
same quantities in the Hartree approximation and found that this simple approximation describes the equilibrium 
properties of the theory remarkably accurately, to within 10% for most observables. We point out that in the high 
temperature, or equivalently the large energy density regime there is a strong rcnormalization of the single particle 
frequencies. 

After studying the equilibrium properties, we presented a method to solve the equations of motion based on the 
dynamics on a light-cone lattice. This method is particularly suitable to study the evolution in field theories within 
the setting of heavy ion collisions where the initial dynamics is mainly along the light cone, it is very stable and 
conserves energy to high accuracy. The equations of motion for the classical field where solved with a broad range 
of initial conditions and energy densities, corresponding to microcanonical evolution. In all of these initial conditions 
the initial energy density was stored in few (or a narrow band) of long wavelength modes. 

Our main results reveal several distinct stages of evolution with different dynamical features: a first transient stage 
during which there are strong fluctuations is dominated by the interaction term, mode mixing begins to transfer 
energy from small to larger wavevectors. A second stage is distinguished by the onset of a very effective transfer 
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FIG. 17: The ultraviolet cascade for k > and three different values of the energy density (top and bottom- left) . In all three 
cases and for all different times, a very good data collapse is obtained, upon rescaling, on the universal scaling function F(u, 0) 
(bottom-right). 



of energy to larger wavevectors, this is an energy cascade towards the ultraviolet. During this second stage 
the interaction term </> 4 diminishes and eventually becomes smaller than the spatio-temporal gradient terms. At the 
end of this stage 4> A <C <fi' 2 ; it 2 . The third stage is dominated by the ultraviolet cascade and the power spectrum 
of the canonical momentum |7Tfe| 2 features a universal scaling form. The cascade is characterized by a front k(t) 
which slowly moves towards the ultraviolet cutoff as k(t) ~ (E/L)* t a with a ~ 0.25 an almost universal exponent 
independent of the lattice spacing and the details of initial conditions but weakly dependent on the energy density. 
During this stage and while the front of the cascade is far away from the cutoff, the power spectrum is universal. 
Behind the front of the cascade the power spectrum is that of thermal equilibrium with an effective temperature 
T e g(t) which slowly decreases towards the equilibrium value. During this stage we find that all observables have 
the same functional form as in thermal equilibrium but with the time dependent effective temperature. Namely the 
wake behind the front of the ultraviolet cascade is a state of local thermodynamic equilibrium. This stage of 
universal cascade ends when the front is near the cutoff, at a time scale t\ ~ -g aT^ with a the lattice spacing. 

The dynamics continues for t > t\ but is no longer universal, true thermalization is actually achieved in the infinite 
time limit and occurs much slower than in any previous stage. 

Thus we find that thermalization is a result of an energy cascade with several distinct dynamical stages. Universality 
and scaling of power spectra emerges in one of the later stages and local thermodynamic equilibrium is associated 
with this stage. While we find scaling and anomalous dynamical exponents we did not find any turbulent cascade. 

We concur with previous studies that thermalization is indeed achieved but on extremely long time scales, however 
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FIG. 18: Numerical results for the universal function F(u, v) obtained through time averages of P(k, t) from t — 51200 to 
t — 384000, for several values of the UV cutoff, at fixed energy density. When a — 0.01 the cascade is almost over and we 
obtain a profile very close to F(u, 1) [see ea . (14. lUt ] ■ When a = 0.005 the bulk of the cascade reaches the UV cutoff, causing 
significant deformations. In the other cases the bulk of the cascade is still far away from the cutoff and F (it, 0) is approached. 
F(u,0) is obtained in the continuum limit a — > 0. 



our study reveals that virialization starts to set earlier than local thermodynamic equilibrium. 

It is important to emphasize that the universal cascade described above is different from a turbulent regime described 
in ref. 29] where the classical massless 4> A model was studied in 3 + 1 dimensions. In that reference, the spectrum for 
the distribution of particles defined with respect to the Hartree frequencies is studied and the numerical results are 
interpreted in terms of wave turbulence and fit to a Kolmogorov spectrum in p9l | . 

We do not find a turbulent spectrum in the power spectra of the field or its canonical momentum, instead we 
find local thermodynamic equilibrium with an effective time dependent temperature which diminishes slowly. While 
our study in 1 + 1 dimensions reveals local thermal equilibrium in contrast to the results of ref. 29] in the 3 + 1 
dimensional case, there are several features in common: a)extremely long thcrmalization time scales, as compared 
to the two natural time scales mT 1 or T _1 , b) a cascade of energy towards the ultraviolet and a cascade front that 
evolves towards the cutoff as a function of time. This latter feature can be gleaned in fig. (2) in ref. pflj . 

The theory of weak wave turbulence which describes the dynamical evolution in terms of kinetic equations leads to 
a turbulent spectrum for nonrelativistic and nonlinear wave equations as the three wave equation 28] . It was argued 
in ref. [30| that the non-equilibrium evolution in cj) 4 in 3 + 1 dimensions but in the small amplitude regime features a 
scaling behavior and a power law albeit different from the one with a ~ 0.21 — 0.25 that we find in 1 + 1 dimensions. 

More recently [3l]], a study of local wave turbulence for a system dominated by four- wave interactions, again via 
wave kinetic equations and following the methods of ref. [2^], reports the formation of cascades that feature a front 
that moves forward in time, the wake behind the front features a Kolmogorov-Zakharov spectrum of weak turbulence. 
The behavior of a moving front found in ref. |3l| is similar to the front of the cascade that we find numerically, but 
the wake is different, we find local thermodynamic equilibrium, while in ref.|3l| the distribution is of the Kolmogorov 
type. 

While a kinetic description may not reliable during the first stage of the dynamics, a suitable kinetic description 
in terms of a distribution function for particles of the (strongly) renormalized frequency may be available during the 
second and third stages. We comment on these and other issues below. 



Discussions and comments: 
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FIG. 19: Numerical estimate of the function h(s) obtained by interpolating the data of a k(t) vs. s = a t 1 ^. Deviations occur 
here when the time t is not large enough compared to the averaging interval and are not included in the interpolation. 

As discussed in section II, the classical limit must be understood with a lattice (or ultraviolet) cutoff to avoid the 
Rayleigh-Jeans divergence. For a finite energy density the naive continuum limit leads to a vanishing temperature 
in equilibrium. Furthermore the classical approximation is valid when occupation numbers are large. In the cases 
under study for large energy density and initial conditions for which the energy is stored in long-wavelength modes, 
the classical approximation is valid and reliable for small wavevectors. Our study reveals that the cascade of energy 
proceeds very slowly after the first stage, thus suggesting that the dynamics after the first stage could be studied 
within a kinetic approach. However, such kinetic approach should include the strong renormalization of the single 
particle frequencies in order to define the slowly varying distribution functions. 

Another result of our study is that the simple Hartree approximation is in remarkably good agreement with the 
exact results in equilibrium. 

Thus we conjecture that the following scenario for studying thermalization in a quantum field theory may prove 
suitable to understanding the dynamics from initial conditions that entail large occupations for long- wavevectors: 
study the numerical evolution in the quantum Hartree approximation up to the time scale at which the interaction 
and consequently the evolution becomes on long time scales. At that stage read out the occupation numbers for 
suitably defined single quasiparticles. Obtain the kinetic equations for the evolution of these quasiparticle occupations 
for a weakly coupled theory, and continue the evolution via these kinetic equations towards final equilibration using 
the occupation numbers at the end of the Hartree evolution. 

This program, if proven suitable, is certainly much more economical and efficient than trying to implement more 
sophisticated and calculational intensive approximation schemes that lead to non-local update equations. This sug- 
gestion is analogous to that advocated in 32] for the evolution after nucleus-nucleus collisions. 

We highlight that we do not find a thermalization threshold as found in the study of FPU chains 0,0. Namely, 
even for low values of the energy density the tp 4 theory reaches thermal equilibrium for very long times. Indeed, the 
thermalization scale t\ diverges in the a — ► limit [see ea. (l4.13|) ]. The fact that the thermalization slows down when 
k(t) approaches the UV cutoff (Brillouin zone) can be understood as follows. The UV cascade can be understood as 
follows. The equation of motion 12.5( 1 for the Fourier component of the classical field with wavevector k features the 
non-linearity ^2 kl X)fc 2 "fei ( t ) k 2 4 > k-k 1 -k 2 - Thus, if initially the amplitude corresponding to a wavevector large k is 
small, the modes with small k and large amplitude lead to a large non-linearity that acts as a source for the amplitude 
with k. In this manner, if the initial power spectrum is localized at small wavevectors the amplitude of the higher 
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FIG. 20: Data collapse over the scaling function h(s) of ak(t) for E/L — 520 and E/L = 104. 
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FIG. 21: The double scaling function h(s) reconstructed from several set of data at different E/L and a. 



modes increases via the non-linearities and the transfer of energy by mode mixing. By periodicity the sum k — k\ — k^ 
is constrained to be in the first Brillouin zone. Therefore when k—k\ — k2 is larger than the cutoff, periodicity implies 
that this sum must be folded back inside the first Brillouin zone by adding a vector in the reciprocal lattice, namely 
7r/2a. Hence once the amplitudes for wavevectors near the cutoff are beginning to be large, the non-linearity transfer 
the energy to low momentum modes. 

Hence, the cascade mechanism is less efficient for k < A than for k <C A. 
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FIG. 22: Effective thermalization of !jP{t). 

While we have focused on the study of the dynamics in a classical field theory with a cutoff as a timely and 
interesting problem all by itself, the applicability of our results to the realm of quantum field theory descriptions of 
early Universe cosmology and ultrarelativistic heavy ion collisions must be discussed. 

Classical field theory is likely to describe the QFT non-equilibrium dynamics at least in the cases in which the 
occupation numbers are very large. Strongly out of equilibrium when a temperature cannot yet be defined, the classical 
description is valid for the range of momenta corresponding to large occupation numbers. Near equilibrium this would 
be the case for frequencies much smaller than the temperature. Our study reveals that effective thermalization is 
reached for times t > to with a time dependent effective temperature. Hence, our classical results should apply to 
QFT for modes with frequencies much smaller than the effective temperature. 

In order to extract a criterion for the validity of the classical approximation we must revert to the physical variables 
T p , a p related to the rescaled variables used in this study as in eqs. 12.2p - Q2.8fl . namely T p — m 3 T/A, a p = 
a/m, k p = mk. T p is the temperature in physical units and the physical momenta in the lattice are restricted to be 
— 7r/a p < k p < 7r y 'a p (first Brillouin zone). 

The exact equilibrium results obtained above indicate that the exact equilibrium frequencies are of the form u)(k) = 
\Jk 2 + bT 2 / 3 with b of order 1, where T is the dimensionless rescaled temperature T — A T p /m 3 . While in the classical 
theory the mass and coupling constant can be absorbed in an overall rescaling of the dimensionful parameters and 
the amplitude of the fields, that is not the case in the quantum theory where the amplitude of the fields and their 
canonical momenta are fixed by the commutation relations. Thus a comparison of the classical and quantum theory 
requires reverting to the physical variables. Moreover, the dimensionless ratio A/m 2 plays a role in the discussion as 
we show below. 

By following the numerical evolution of single modes it is easy to see that for all modes in the bulk of the cascade, 
where local equilibrium holds, the time-dependent dimensionless frequencies have the expected form 

- V / fc 2 + 6T cff (0 2 / 3 

We emphasize that the cascade corresponds to a state of local equilibrium with the effective temperature T e f f (t) which 
is much larger than the final thermodynamic equilibrium temperature. 

A simple criterion for the validity of the classical approximation in QFT is that 



Tp,eK(t) > LU p (kp(t)) 



(5.1) 
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FIG. 23: Effective thermalization of 4 (i). 



where uj p (k p ) — \/k p + b [A T PjC ff (i)] 2 / 3 are the dimensionful frequencies. In dimensionless variables ea. (|5.1f) takes the 
form, 



rn 



Teff(t) » m Jk{t) 2 + &T eff (i) 2 / 3 



(5.2) 



It is enough to discuss the modes at the wake to see the applicability of the classical approximation since k < k(t) for the 
modes in the cascade. At times early enough, in regimes such that the plateau is already formed but k{t) <C T c s (t) 1 / 3 , 
ea. (|5.1|) holds provided T e fi{t) ^> (A/m 2 ) 3 / 2 , which is certainly true for A small enough. 
On the other hand if k(t) ^> Teff(i) 1 / 3 the classical approximation will be valid provided 



m 

T 



Teff(i) » k{t) 



Using our late time estimates eqs. (14.811 . (I4.13|) for a = 0.25 one obtains the condition 

m 4 E En 



t <C tQ , tQ 



A 2 L XL r 



(5.3) 



(5.4) 



This shows that for sufficiently small coupling constant and fixed UV cutoff A = 7r/a p the classical approximation 
could be valid up to full thermalization, that is tQ ^> t%, which using the first of eq. (|4.13Jl is equivalent to the 
condition 



E A 

> 



Recalling that E/L ~ T/(2a), this is nothing but T p a p ^> 1. In fact, as the front of the cascade moves towards the 
cutoff, for a fixed cutoff, the classical approximation will continue to be valid for all modes behind the wake of the 



cascade if T p a p 3> 1 . 



The validity of the classical approximation for T p a p 3> 1 is well known in the elementary theory of harmonic solids 
where phonons are the elementary excitations. In the Debye model of solids the classical limit is valid for T p > Tu 
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FIG. 24: Effective thermalization of \<f>\ 2 (k,t) for different values of k <C 2k(t) . 



where Tp oc 7r/g„ is the Debye Temperature (and tt /a is the edge of the first Brillouin zone, namely the maximum 
momentum) • 

As per the discussion in section III II the non-linearities are subdominant for Ta 3 <C 1, condition which all but 
guarantees the cascade [see discussion in section lrTT] . Therefore, the validity of the classical approximation is warranted 
provided, 

T„a„ 1 or equivalently Ta^$> — » , (5-5) 

ni 

which is the same condition for the validity of the classical limit in elastic solids. 

Since our results are valid for Ta 3 < 1 the classical results obtained here apply in the quantum theory for a fixed 
cutoff in the whole dynamical range where eq. l|5.3[) is valid. 

We have studied the dynamical evolution in a broad range of parameters from 0.001 < Ta < 1 (with Ta 3 <C 1) and 
found a similar behavior for the dynamics during the stage of the universal cascade. Therefore, our classical results 
applies in QFT provided one chooses values of A/m 2 fulfilling ea. l|5.5|l for the values of Ta indicated in the respective 
plots. Indeed, our classical results turn to be relevant for the weakly coupled quantum field theory. Otherwise, quantum 
effects certainly appear in the process of thermalization [TEl Hsj . Nonetheless, even in the continuum limit at fixed 
coupling constant, when inevitably T p a p — > 0, it is conceivable that, as long as t <C £q, the classical approximation 
is quite accurate for all modes in the bulk of the cascade. One should expect that the shape of the cascade at the 
forefront indeed depends sensibly on quantum effects, since occupation numbers fall off rapidly across it. But the rate 
at which the plateau behind decreases is most likely almost purely classical, at least until occupation numbers are 
high, that is for t <C £q, when the effective temperature is high enough. 

Hence, in the quantum field theory with initial conditions strongly out of equilibrium with a power spectrum 
featuring large occupations for long wavelengths we expect that the initial stages of the dynamics will likely be well 
described by the classical field theory. In this case we expect a cascade to form, the modes behind the front will be 
well described by the classical field theory, whereas in the theory with a large cutoff, the modes with small occupation 
in front of the cascade all the way to the cutoff will require a quantum description. The backreaction of these quantum 
modes onto the long-wavelength modes will modify the time evolution of the cascade. 

In the quantum 4 theory we expect that the time scale for complete thermalization to be finite even in the 
continuum limit. While in the classical limit the equilibrium power spectrum for the field falls off with a power law 
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FIG. 25: The ratio ^rr(i) as a function of the logarithm of the time t for E/L = 3424, 50.28 and 3.669 for L = 20 and 
a — 0.01. This ratio in thermal equilibrium takes values between the low temperature limit 1/3 and the high temperature limit 
0.37077 . . . depicted as a dotted line. Hence, effective thermalization may start only when the ratio falls within such strip [see 
cq . 114. 16^ 1 . For E/L = 3.669 effective thermalization only starts for hit > 20. 

and that of the canonical momentum is flat, hence sensitive to the cutoff, in the quantum theory the exponential 
suppression will make the theory insensitive to the cutoff in the limit Ta<l. Hence, we conjecture that time scale 
for complete thermalization diverges in the continuum quantum theory when H — > 0. Furthermore, we expect faster 
thermalization in the quantum theory since zero point fluctuations imply finite amplitude for all modes in the initial 
state. 

The numerical studies previously reported on the dynamics of classical and quantum field theories |ll | -|20 | have 
not yet focused on studying the mechanism of energy transfer from long to short wavelengths as a function of time, 
although ref.[2(j reported the emergence of spatio-temporal structures. It would indeed be interesting to study if 
the universal cascade found in the classical theory remains at least during some early and intermediate stages in the 
quantum theory. 

Our study in conjunction with a previous body of work clearly suggests that thermalization, at least in this 
scalar theory, occurs on extremely long time scales, much longer than the natural scales m _1 , T^ 1 . However, local 
thermodynamic equilibrium with a slowly decreasing effective temperature T D g (t) emerges relatively soon. This result 
all by itself raises important questions on the current assumptions on thermalization both in cosmology as well as in 
heavy ion collisions. If thermalization involves a time dependent effective temperature, both in cosmology as well as 
in the formation and evolution of the quark gluon plasma, there are immediate phenomenological consequences and 
a reassessment of the current assumptions and models is called for. 

An intense program to study these issues and their impact on both fields is therefore warranted. 

VI. APPENDIX 

We provide here more details on our numerical calculations. They were performed on a cluster of personal computers 
at LPTHE-Paris and at the Physics Department of Milano-Bicocca. The computers were used in parallel with different 
values of relevant parameters on each one of them. 

In our numerical evolution the number of iterations of the discrete system, ea. (|3.6|l . is t max /(2a) when evolving 
from t = up to t — t max - The computer time increases for smaller spacing also because at fixed size L, smaller a 
implies larger N — L/ {2a). This implies a number of operations of order a~ 2 . Moreover, if fast Fourier transforms are 
performed every iteration, to average the power spectra over the time intervals [t, t + r] covering completely the time 
span from t = up to t = i max , then altogether the computer time scales as a -2 log - for small a. Actually, when 




44 



also averages over initial conditions are performed, it turns out that r need not be so large as to provide a complete 
covering, since the average over initial conditions is much more effective in reducing fluctuations than averaging in 
time [see fig ES]- In other words, it is sufficient to consider time intervals [t,t + r] that cover only a relatively small 
part of the full time span from t = until t — t max . This implies that computer time really scales just as a~ 2 for 
small a. 

As pointed out in sec. Iff! El we have two methods to compute \<fik(t)\ 2 - Either, to extract the field <fi(x,t) from the 
lattice fields F(n, s) and G(n, s), Fourier-transform it to (j>k(t) and then perform all needed averages on \<j)k{t)\ 2 . Or 
to directly compute the averages of the correlations of F(n, s) and G(n, s) and extract from them the correlations 
4>(f>(x, t) and 7f7f(x, £). Thanks to the exceptional scaling properties of the Fast Fourier Transform, the first method is 
actually preferable to the second, if the space average in ea. (|3.f4|l is really performed and very extensive time averages 
are considered. However, as just explained, the time averages can be greatly reduced when averaging on numerous 
enough initial conditions. 

We have seen in the previous sections that the mean wavenumber k(t) plays a central role in the description of the 
UV cascade. However, the numerical determination of the exponent a and of the amplitude hi for the time behaviour 
of k(t) in the scaling window, ea. l|4.3|) or, even more generally, in the determination of scaling function h(s) in ea. (|4.3|) 
is a nontrivial matter. This is due to the characteristic fluctuations of the data that limit the control of the sistematic 
effects caused by the initial transient, by the scaling violations due to the UV cutoff and by the averaging over time 
intervals when the time is not much larger than interval length r. This type of fluctuations can be systematically 
reduced by enlarging the set of initial conditions over which we average, but are still present in a significant way for 
the number (typically 20 to 30) of initial configurations we used. 

We identify the initial transient in k(t) by a self-consistent fitting procedure. Suppose a very good fit has been 
obtained for the data points at times larger than a certain time t (the actual fitting procedure is described below). 
Then we can see if and how much the data at t < t deviate from the fitting function. If the fit starts to worsen 
immediately before t and in a steady manner towards t — 0, then the transient cannot be fully disentangled from the 
fitting errors and t needs to be increased until we find that the fit remains quite good for a while before t and then 
worsen in rather abrupt way. The time at which this worsening occurs provides an alternative estimate of the time 
scale to introduced in sec. IIIIDI to signal the end of the transient and the start of the UV cascade that will become 
universal for t 3> to. Using the parametrization 

k(t) = tnt a [1 + <?(*)] , 

can then (approximatively) read to out of the plot of g(t) vs. t [see fig 12 7|. We see that log to ^ H (in agreement 
with what found in sec. IIII Dl for other observables) more or less independently on E/L as long as E/L > 52. 
When E/L = 2.6, to is much larger instead, logio ~ 15, as found also for the other observables. Notice that, when 
E/L = 2.6, k(t) is very well fitted for t < t < t max , with log? ~ 14.5 and by logi max ~ 17.5 by the pure power law 
1|4.3(1 of the scaling window, but with a power a ~ 0.186 definitely smaller than those at higher energy densities. 

Having disposed of the initial transient, we now turn to the fitting procedure to determine h\ and a. As a matter 
of fact, the simple, direct fitting of the numerical data of average wavenumber k(t) with ea. H4.3|) is not quite precise, 
even when the scaling window is identified. For instance, the almost linear behaviour in the log-log scale for smaller 
values of the spacing a is evident from fig. 1161 but the actual values of a and of hi are poorly fixed in this simple 
fitting approach. 

To give an idea of the situation, let us consider a specific example. In a parallel run of 20 evolutions from t — to 
t = 1.28 . . . 10 7 , with a = 0.0025 and E/L = 520, each one starting from the same initial wavenumber k^s filled, but 
with randomly chosen amplitude and phases [according to the discussion in sec. IHIB] . a 20 x 400 matrix M of data 
for k(t) is obtained. Each row of M corresponds to a given choice of initial conditions and the n— th column contains 
the values of the mean wavenumber at t — n * At for n — 1,2, ... , 400 and At = 32000. In this runs the time average 
intervals are rather small (r = 100) compared to the data taking interval At and so they cause a negligible sistematic 
error. On the other hand the effects of the finite UV cutoff are not so small [logfc(t) bends downward as a function 
of log t] and this facts suggests to try a simple least squares fit with 

Te(t) ~ hi t a + h 2 a t 2a , (6.1) 

using hi, /12 and a as free fitting parameters. When the fit is performed on the mean wavenumber extracted from the 
power spectrum |7r| 2 (fc,i) averaged over the initial conditions, we obtain 

hi = 4.3942 , h 2 = -7.4101 , a = 0.2401 . (6.2) 

This differs very slightly from the result of the fit over the column-averaged M, that is 



hi = 4.4200 , h 2 = -7.4734 , a = 0.2397 



(6.3) 
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TABLE I: Numerical fits for the expansion coefficients of the scaling function h(s). Relative statistical errors on the data are 
of the size of few percent on hi , but unavoidably larger for or hs , since these correspond to corrections that are too small 
in most cases. When /12 or /13 are not included in the fit, an asterisk is shown. Another source of (sistematic) errors is the 
approximate validity of the scaling hypothesis ea. 14.31 for a not small enough. 

since k(t) depends in a weak nonlinear way on |-7r| 2 (fc,t). However, to give an estimate of the statistical error on this 
figures, we should repeat the procedure on each row of M separately. This yields 

hi =4.7660 ±1.7855, h 2 = -1.9006 ± 11.3948 , a = 0.2373 ± 0.0337 . (6.4) 

Of course the mean values here do not coincide with the values in either eq. (|6.2() or eq. (|6.3|l since the fitting parameters 
depend very nonlinear ly on the data. The annoying point however is the large statistical error on this figures, which 
significantly exceeds the differences in the mean values. In principle, relying on the standard error reductions in the 
mean of independent random variables, we could assign error bars ea. il6.2|) or ea. l|6.3fl smaller by a factor 1/V20 with 
respect to those in ea. l|6.4[) . However, at the same time the confidence level of each separate fit on the rows of M is 
much lower than in that yielding the estimates in ea. l|6.2|l or ea. l|6.3[l . since the larger fluctuations (more or less by 
the inverse factor y/20) make much larger the least squares of the fit. Altogether this entail a relative error ranging 
from 10% for h\ and 5% for a (with the most optimistic attitude) to 40% for hi and 20% for a in a more conservative 
approach. The situations for h 2 is much worse, but it does not have the same relevance of hi and a in the continuum 
theory. 

Interestingly enough, it turns out that the full application of scaling hypothesis of ea. l|4.3|l is more efficient that 
direct fitting in pinpointing the value of hi and especially a. That is, the numerical data of a k(t), for different values 
of a at fixed E/L, collapse indeed very neatly on a single curve when plotted versus s = a t a (thus verifying the 
scaling idea), if a has a specific value. Small changes of a, of relative order 10~ 2 , spoil the collapse in a evident and 
easily quantifiable fashion. In this way we determined with few percent accuracy the values a — 0.222 for E/L = 104, 
a = 0.231 for E/L = 520 and a = 0.25 for E/L = 5200 [see figs. IT51 and I2D] . Next, with a fixed, we can repeat the 
fit with the expansion of a k(t) in powers of a t a as in eas. l|6.1|l . obtaining a much better determination for hi and 
often also for h 2 [see table [Q. In particular, when the fit with fixed a = 0.231 is considered in the example discussed 
above for each row of the matrix M separately, we obtain 

hi = 4.98167 ± 0.1389 , h 2 = -8.3544 ± 1.5271 , 

to be compared with the fit in ea. (|6.4|) when a is a fitting parameter. 

We see that the parameters hi and a do depend on the energy density E/L = 520. a grows from a ~ 0.186 when 
E/L = 2.6 to a ~ 0.25 when E/L = 5200. We did not attempt an analytic fit of a(E/L), but it does seem to saturate 
fast at higher energy densities. 



Acknowledgements: D. B. thanks the N.S.F. for partial support through grants PHY-0242134, and the hospitality 
of LPTHE where part of this work was carried out. C. D. thanks the LPTHE for the warm hospitality. H. J. d. V. 



4G 



3.5 



2.5 



1.5 



0.5 



E/L=104, t = 3200, n = 1 
E/L= 104, x = 3200, n = 20 
E/L = 520, T = 100, n = 1 
E/L = 520, x = 100, n = 20 




8 9 10 11 12 13 14 15 16 17 

iog(t) 

FIG. 26: Comparison of different averagin methods, n is the number of initial conditions in each average. 



thanks the Department of Physics and Astronomy at the University of Pittsburgh for their hospitality. 



[1] E. Fermi, J. Pasta and S. Ulam, Note e Memorie vol. II (Collected Papers) p. 978, Accademia dei Lincei, 1965. 

[2] R. Livi, M. Pettini, S. Ruffo, M. Sparpaglione, A. Vulpiani, Phys. Rev. A28, 3544 (1983), A31, 1039 (1985). M. Pettini, 
M. Landolfi, Phys. Rev. A41, 768 (1990). J. DeLuca, A. J. Lichtenberg, S. Ruffo, Phys. Rev. E51, 2877 (1995), E54, 2329 
(1996), E60, 3781 (1999). J. DeLuca, A. J. Lichtenberg, M. A. Lieberman, Chaos, 5, 283 (1995). 

[3] G. Parisi, Europhys. Lett. 40, 357 (1997). 

[4] E. W. Kolb and M. S. Turner, The Early Universe, Addison- Wesley Pub. Co. Reading, MA, 1990. 
[5] A. R. Liddle and D. H. Lyth, 'Cosmological Inflation and Large Scale Structure', Cambridge Univ. Press, 2000. 
[6] A. H. Mueller, Nucl. Phys. A702, 65c (2002), Nucl. Phys. B572, 227 (2000) and Phys. Lett. B475, 220 (2000). 
[7] R. Baier, A. H. Mueller, D. Schiff, D. T. Son, Phys. Lett. B502, 51 (2001). 

[8] A. Krasnitz, R. Venugopalan, Nucl. Phys. B557,237 (1999), Phys. Rev. Lett. 84, 4309 (2000) and Phys. Rev. Lett. 86, 
1717 (2001). J. Bjoraker, R. Venugopalan, Phys. Rev. C63, 024609 (2001); A. Krasnitz, Y. Nara, R. Venugopalan, Phys. 
Rev. Lett. 87, 192302 (2001), Nucl. Phys. A717, 268 (2003), hep-ph/0305112 
[9] J. Traschen and R. Brandenberger, Phys Rev D42, 2491 (1990), Y. Shtanov, J. Traschen and R. Brandenberger, Phys. 
Rev. D51, 5438 (1995). L. Kofman, A. Linde and A. A. Starobinsky, Phys. Rev. D56, 3258 (1997), G. Felder, L. Kofman, 
A. Linde, Phys. Rev. D59, 123523 (1999); G. Felder, L. Kofman, A. Linde, I. Tkachev, JHEP 0008, 010 (2000). 

[10] D. Boyanovsky, H. J. de Vega, R. Holman, J. F. J. Salgado, Phys. Rev. D54,7570 (1996); D. Boyanovsky, H. J. de Vega, 
R. Holman, D.S-Lee, A. Singh, Phys. Rev. D51, 4419 (1995); D. Boyanovsky, D. Cormier, H. J. de Vega, R. Holman, 
Phys. Rev. D55, 3373 (1997); D. Boyanovsky, D. Cormier, H. J. de Vega, R. Holman, A. Singh, M. Srednicki, Phys. Rev. 
D56, 1939 (1997); D. Boyanovsky, D. Cormier, H. J. de Vega, R. Holman, S. P. Kumar, Phys. Rev. D57 , 2166 (1998); 
D. Boyanovsky, C. Destri, H.J. de Vega, R. Holman, J. F. J. Salga do, Phys. Rev. D5 7, 7388 (1998). 

[11] Sz. Borsanyi, A. Patkos, D. Sexty, Phys. Rev. D66, 025014 (2002), hep-ph/0303147 and hep-ph/0301117 

[12] G. Aarts, J. Smit, Nucl. Phys. B511, 451 (1998), Nucl. Phys. B555, 355 (1999) and Phys. Rev. D61, 025002 (2000). 

[13] G Aarts, G. F. Bonini, C. Wetterich, Nucl. Phys. B587, 403 (2000) and C. Wetterich, Phys. Rev. D63, 025012 (2001). 

[14] G. Aarts, J. Berges, Phys. Rev.Lett. 88, 041603 (2002); G. Aarts, Phys. Lett. B518, 315 (2001). 

[15] J. Berges, Nucl. Phys. A702, 351 (2002) and Nucl. Phys. A699 , 847 (2002). J. Berges, J. Cox , Phys. Lett. B517, 369 

(2001); G. Aarts, D. Ahrensmeier, R. Baier, J. Berges, J. Serreau, Phys. Rev. D66, 045008 (2002). 
[16] K. B. Blagoev, F. Cooper, J. F. Dawson, B. Mihaila, Phys. Rev. D64 125003 (2001); F. Cooper, A. Khare, H. Rose, Phys. 



47 



0.35 



0.25 - 



a = 0.001 25, E/L = 5200 

a = 0.001 25, E/L = 104 

a = 0.001 25, E/L = 520 

a = 0.005, E/L = 52 

a = 0.01, E/L = 2. 



*s 0.15 




0.05 



-0.05 



FIG. 27: Numerical determinations of the function g(t) that describes the initial transient for several values of E/L. 



Lett. B515, 463 (2001); F. Cooper, J. F. Dawson, B. Mihaila, Phys. Rev. D67, 056003 (2003). 
[17] J. Baacke, A. Heinen, hep-ph/0212312 and hep-ph/0305220 

[18] J. Berges and J. Serreau, |hep-ph/03022l"o| J . Berges, Sz. Borsanyi, J. S erreau, Nucl.Phys. B660, 51 (2003); J. Serreau, 

Nucl.Phys. A715, 805(2003); J. Ber ges and J. Serre au, hep-ph/0208070] 
[19] D. J. Bedingham and H. F. Jones, |hep-ph/02090"60"] 
[20] M. Gleiser, R. C. Howell, hep-ph/0209176 

[21] G. Amelino-Camelia and S.-Y. Pi, Phys.Rev. D47,2356 (1993); M. Salle, J. Smit, J. C. Vink , Nucl.Phys. B625495, (2002); 
Phys.Rev. D64, 025016 (2001); L. M. A. Bettencourt, Karen Pao and J.G. Sanderson, Phys.Rev. D65 025015 (2002), and 
references therein. 

[22] R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, Mc-Graw Hill, 1965. 
[23] K. Symanzik, unpublished, circa 1969. 

[24] K. Banerjee, S. P. Bhatnagar, V. Choudhry and S. S. Kanwal, Proc. R. Soc. Lond. A 360, 575 (1978). 
[25] F. T. Hioe, D. MacMillen and E. W. Montroll, Phys. Rep. 43, 305 (1978). 

[26] See for example, C. Destri, H. J. de Vega, Nucl. Phys. B 290, 363 (1987) and references therein. 
[27] F. Zanlungo, Thesis for the Italian 'Laurea in Fisica', Universita Milano-Bicocca, unpublished. 

[28] V. E. Zakharov, V. S. L'vov and G. Falkovich, Kolmogorov Spectra of Turbulence I: wave turbulence, Springer Verlag, Berlin, 
1992; V. S. L'vov, Wave Turbulence Under Parametric Excitation: Application to Magnets, Springer Verlag, Berlin, 1994. 
[29] R. Micha and I. I. Tkachev, Phys. Rev. Lett. 90, 121301 (2003) and hep-ph/0301249 
[30] D. T. Son, Phys. Rev. D54, 3745 (1996). 

[31] C. Connaughton, A. C. Ne well and Y. Pome au, [physics/0304076 
[32] A.H.Mueller and D.T.Son, hep-ph/0212198 

[33] See for example the textbook: R.K. Pathria, Statistical Mechanics, Pergamon Press, N.Y. 1972. 



